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Introduction  and  Background 

^ — - ^  The  DIODE  computer  prograr/at  the  Harry  Diamond  Laboratories 

(HDL)  is  used  to  compute  current  and  voltage  characteristics  in 
semiconductors  or  gases,  has  been  used  without  change  for 
approximately  12  years;  improvements  upon  the  program  are  over¬ 
due.  This  report  documents  the  existing  program  and  will  serve  as  a 
basis  upon  which  the  improvements  will  be  built.  A  listing  of  the 
computer  code  is  provided  in  appendix  A. 

HDL's  interest  in  the  basic  physical  processes  of  electrical  breakdown 
in  gases  arose  from  the  need  for  sensitive  trigger  tubes  for  use  in 
electrostatic  fuzes.  The  initial  theoretical  work  was  based  on  papers  by 
Varney  et  al  [1]  and  by  Crowe  et  al  [2].  In  the  mid-1950's,  Professor 
William  Dow,  Chairman  of  the  Electrical  Engineering  Department  at 
the  University  of  Michigan  and  member  of  the  Scientific  Advisory 
Board  of  the  Diamond  Ordnance  Fuze  Laboratories  (the  predecessor 
of  HDL),  suggested  that  the  problem  could  best  be  solved  on  the  then- 
new  electronic  computers.  Ward  modified  the  formulation  of  Crowe 
et  al  to  apply  to  the  rare  gases  and  used  the  National  Bureau  of 
Standards  (NBS)*  SE  AC  (South  East  Automatic  Calculator)  computer 
to  calculate  results.  In  order  to  extend  the  calculations  to  higher 
current  densities,  Ward  included  the  space  charge  of  electrons  and 
reformulated  and  encoded  the  problem  for  tine  IBM  704  computer. 

Time  dependence  was  included  for  the  first  time  when  Borsch-  Supan 
and  Oser  [3]  wrote  the  basic  program  that  is  still  in  use.  Thanks  to  the 
insistence  of  Irene  Stugun,  then  chief  of  the  NBS  programming  sec¬ 
tion,  the  basic  continuity  equations  were  used,  rather  than  the  partially 
integrated  equations  used  previously.  Over  100  papers,  reports,  and 
oral  presentations  (see  bibliography)  have  resulted  from  the  DIODE 
program. 

In  the  summer  of  1965,  a  summer  student  at  HDL,  Edward  R.  Berman, 
modified  the  program  to  include  semiconductors  [4].  The  major 
changes  were  the  inclusion  of  ionization  by  holes  and  the  use  of  double 
precision  because  of  the  higher  densities  present  in  semiconductors. 
In  1969,  under  the  direction  of  Burton  Udelson,  Howard  Bloom  [5] 
modified  the  external  circuit  of  DIODE  in  order  to  study  IMP  ATT 
oscillators.  Bloom's  report  includes  a  full  program  listing.  Since  that 
time.  Ward  [6]  added  thermal  effects  to  the  program.  Until  1988, 
computer  cards  were  used  for  the  input  data.  Now  the  program  may 
be  used  in  the  time-sharing  mode  at  a  remote  terminal. 

Before  the  computer  code  and  its  use  are  discussed,  we  present  the 
mathematics  that  forms  the  basis  of  the  program.  We  then  discuss  in 
detail  the  format  of  the  input  data  required  to  run  the  program;  the 

*Now  the  National  Institute  for  Standards  and  Technology  (NIST). 
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input  data  required  consist  of  physical  properties  of  the  material 
studied  and  device  descriptions.  The  information  must  be  supplied  by 
the  user;  HDL  users  can  build  on  datasets  previously  used,  changing 
only  the  information  necessary.  The  output  produced  by  the  program 
is  then  described  in  detail  (an  example  of  the  output  is  provided  in 
app  B). 

As  an  aid  to  users  in  devising  input  datasets,  we  include  information 
on  the  material  parameters  that  have  been  commonly  used;  tins 
information  includes  both  data  and  suggestions  on  how  to  choose 
constants  cased  on  various  earlier  efforts. 

Finally,  we  discuss  the  interactive  use  of  DIODE  as  it  has  been 
implemented  on  HDL's  IBM  3090  mainframe.  These  instructions  are 
particular  to  HDL  users,  but  non-HDL  users  may  find  them  useful  if 
they  intend  to  implement  the  program  elsewhere. 

2.  Formulation  of  Computer  Program 

2.1  Basic  Equations 

The  one-dimensional  continuity  equations  for  electrons  and  holes, 
respectively,  in  a  semiconductor  are 

^  =  ^  +  ,  (1) 

at  ax 

^=^+c^.+®.-R  ,  (2) 

where  a:  and  f  are  the  space  and  time  coordinates,  e  is  the  electron  charge, 
n  and  p  are  the  electron  and  hole  number  densities,  a  and  P  are  the 
ionization  coefficients  for  electrons  and  holes,  and  R  is  the  recombi¬ 
nation  rate.  The  electron  and  hole  current  densities, /_  and  /^,  are  given 
by 

l  =  neiii.E-eD.[^]  ,  (3) 

U  =  pefuE-eD^i^^]  ,  (4) 

where  IJ._  and  are  the  electron  and  hole  mobilities,  D_  and  are  the 

diffusion  coefficients,  and  £  is  the  electric  field. 

Space-charge  effects  are  determined  by  the  one-dimensional  Poisson 
equation 

dE_e{n-p+N)  ^  ^5. 

dx  ^ 
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where  Kis  the  dielectric  constant  for  the  semiconductor,  is  the  per¬ 
mittivity  of  free  space,  and  N(x)  is  the  distribution  of  net  fixed  charge. 
The  sign  convention  in  equations  (1)  through  (5)  is  chosen  so  that  E, 
/_,  and  are  all  normally  positive  quantities. 

The  one-dimensional  thermal  diffusivity  equation  is 


I. 


(6) 


where  T  is  the  temperature,  /  is  the  total  current  density,  p  is  the 
semiconductor  density,  c  is  the  heat  capacity,  and  kj  is  the  thermal 
conductivity. 

It  may  be  shown  from  equations  (1)  through  (5)  that  the  total  current 
density, 

(7) 

is  a  constant  in  space.  This  merely  expresses  current  continuity  in  one 
dimension.  The  accuracy  of  the  calculations  may  be  monitored  by  the 
constancy  of  /  across  the  width  of  the  diode. 

The  variation  of  the  intrinsic  density  with  temperature  was  chosen  as 


«?(!)= 2  X  exp||j  (SOOr 


-r^) 


(8) 


where  Eg  is  the  bandgap  energy  of  silicon.  The  variation  of  the  injected 
(thermally  generated)  current  density  with  temperature  was  assumed 
to  be  the  same  as  for  n, ^7),  since  np  =  nf  and  the  majority  carrier  den¬ 
sity  is  fixed  by  the  doping  level. 


2.2  Supplementaiy  Equations 

The  cathode  (x  =  0)  to  anode  {x  =  d)  distance  is  divided  into  M  equal 
intervals  of  width  Ax.  Initial  (t = 0)  arrays  of  M + 1  values  must  be  given 
for  n,  p,  and  T,  and  a  similar  array  given  for  the  fixed  charges,  N  = 

-  where  is  the  net  number  density  of  donors  and  Ny^  is  the  net 

number  density  of  acceptors. 

Two  boundary  conditions  are  required  for  equations  (1)  and  (2).  the 
electron  current  density  at  x  =  0  and  the  hole  current  density  atx  =  d. 
Alternatively,  the  number  densities  can  be  given  and  the  current 
densities  calculated  from  equations  (3)  and  (4).  At  present,  the  bound¬ 
ary  currents  are  constant  in  time. 

The  boundary  condition  for  equation  (5)  is  supplied  by  the  total 
voltage  across  the  diode.  The  initial  voltage  across  the  diode  must  be 
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given,  and  the  voltage  for  later  times  is  determined  by  the  external 
circuit.  Optional  external  circuits  are  available,  but  for  this  report,  the 
diode,  shunted  by  a  capacitance,  C,  is  in  series  with  a  load  resistance, 
Rg,  and  a  voltage  source,  V{t).  The  voltage  source  may  be  constant  in 
time,  have  one  discrete  step,  or  have  an  incremental  sinusoidal  varia¬ 
tion.  The  last  option  allows  a  constant  dV/dt  value  to  be  closely 
approximated. 

3.  Input  Data 

The  input  data  for  DIODE  are  formulated  on  21  lines,  preceded  by  a 
title  line,  and  followed  by  a  variable  number  of  doping  distribution 
and  possibly  temperatme  distribution  lines.  The  final  line  is  an  ending 
indicator.  Tlie  data  field  is  justified  right  and  uses  standard  Fortran 
criteria  for  digital  data  and  integers.  The  standard  line  uses  five 
number  fields  of  14  spaces  each;  i.e.,  numbers  end  on  spaces  14, 28, 42, 
56,  and  70.  Care  must  be  taken  that  exponents  end  at  the  correct 
position,  since  zeros  are  read  in  blank  spaces.  Some  lines  have  integers 
both  in  space  1  and  in  space  80. 

Table  1  lists  the  21  lines  with  a  short  identifying  name.  In  the  field 
width  column,  an  integer  in  space  1  is  indicated  by  "I"  and  in  space  80 
by  "I*."  In  the  descriptions  that  follow,  consult  table  1  to  see  the  format 
and  the  array  of  parameters. 

Title  line.  On  the  title  line,  the  first  character  must  be  in  space  1;  a 
maximum  of  71  characters  is  allowed.  Typically,  such  information  as 
the  date  of  the  run  and  the  diode  characteristics  would  be  entered  in 
the  title  line. 

Lines  1  to  3.  On  lines  1, 2,  and  3,  enter  the  avalanche  (that  is,  ionization 
by  collision)  coefficients  for  electrons,  a  (parameters  A1  to  A5);  holes, 
P  (parameters  B1  to  B5);  and  gas  excitation,  d  (parameters  D1  to  D5); 
respectively.  (Since  the  format  is  the  same  for  each  of  these  lines,  only 
line  1  is  described.)  When  a=0,  alpha  parameter  MODA  should  be  set 
at  zero,  and  A1  and  A5  are  irrelevant.  MODA  should  be  set  to  1  if 

a=pA\  exp  /  if  ®  <  A5  ,  (9) 

a=p  A3  exp  /  if  ®  >  A5  ,  (10) 

where  p,  the  gas  pressure,  is  set  equal  to  1  for  semiconductors,  and  E 
is  the  electric  field. 

If  rare  gases  are  being  used,  MODA  should  be  set  to  2;  in  this  case, 
(l£|/p)^/2  replaces  |£|/p  in  equations  (9)  and  (10), 


8 


f2 


(/> 

•a 

bi 

Si 

a 

(X 

B 

o 

U 


VO 


ID 


CD 


cs 


<U  01 
2  60 
I  S 

!3;S 
tx  ^ 


®  « 
■S  2 

TS  OJ 


<ii 

6 

« 

2 


<u 

.5 

■J 


1  I 


I  I 


I 

1  < 


CQ 


I  ^  2 


I  ^ 

I  < 


< 

D 

O 


m 


CD 

8 


Tj<  Tl( 


(Q 

0)  X 

•53  ^ 

P  < 


r-l  C-J 


9 


means  that  first  number  infield  is  an  integer, 
^l*  means  that  number  in  column  80  is  an  integer. 


The  program  computes  A6  internally  to  ensure  continuity  in  a  at  A5; 
the  result  is  included  in  the  printed  output. 

Since 
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A6  =  -A5  [In  A1  -  In  A3  -  A2/A5]  , 
one  must  not  set  A5  equal  to  zero.  MODA  should  be  set  to  3  if 

a  =  pAl  exp  j=^-j  +p  A3  exp  . 

If  MODA  is  set  to  4,  (p/\E\)^^^  replaces  p/\E\. 

Lines  2  and  3,  for  holes  and  excitations,  respectively,  pertain  to  the 
same  equations  as  line  1  (for  electrons). 

Line  4.  On  line  4,  you  supply  data  needed  for  temperature  effects.  If  the 
temperature  is  constant  in  space  and  time,  the  temperature  parameter 
MODI  should  be  set  to  0;  put  the  temperature  value  in  TEMPK,  the 
last  field  in  the  line.  If  you  calculate  that  the  temperature  will  change 
as  a  result  of  power  dissipated,  set  MOOT  =  1.  In  CAT,  put  the 
fractional  change  in  A1  for  a  one-degree  increase  in  temperature,  and 
in  CBT  put  the  fractional  change  in  Bl.  In  PWRM,  put  the  power-law 
dependence  of  the  low  field  mobility  of  electrons  upon  temperature, 
i.e.,  tnJX)  =  tK^(300)T‘’'^^.  Likewise,  in  PWRP,  put  the  exponent  for 
the  hole  mobility.  If  you  are  generating  an  initial  distribution  of 
temperature  after  the  doping  concentration  lines,  MODT  must  be  set 
to  2  (see  discussion  of  "Other  lines"  at  the  end  of  this  section). 

Line  5.  On  line  5,  you  enter  the  material  parameters  of  density,  specific 
heat,  and  thermal  conductivity.  These  parameters  do  not  change  with 
temperature  in  the  DIODE  program. 

Lines  6  to  10.  On  lines  6  through  10,  you  enter  data  describing  the 
variation  of  the  electron  and  hole  mobility  with  field  and  doping  level. 
If  the  hole  mobility,  Pp,  should  be  calculated  according  to 

/^  =  MUP(1-C1£)  ,  £<C3  (11) 

Pp  =  C2  £-^/2  (1  _  C9  £)-3/2  ^  £  >  C3  ,  (12) 

MODVP  should  be  set  to  0  (line  6).  The  computer  calculates  C9  to 
ensure  continuity  of  Pp  at  £  =  C3,  and  C9  is  printed  out. 

If  you  set  MODVP  to  1,  you  must  enter  C9  and  CIS  in  line  10.  Then  the 
computer  solves  equations  (11)  and  (12)  to  obtain  continuity  in  Pp  and 
dpp/dE  over  the  range  from  £  =  C3  to  £  =  CIS  for  the  power  law 

Pp  =  cn+  C12£  +  C13-£2  +  C14£3 

and  prints  out  Cll  to  C14. 


(13) 


The  last  parameter  on  line  6,  C4,  is  the  maximum  hole  velocity,  usually 
termed  saturation  velocity,  to  be  used  in  the  field  range  for  equation 
(12).  Since  a  peak  velocity  at  low  fields  can  exceed  the  saturation 
velocity,  as  is  the  case  for  electrons  in  GaAs,  the  maximum  hole 
velocity  does  not  apply  for  fields  less  than  C3  (eq  (11)). 

If  you  set  MODVNP  to  0,  in  line  7,  the  hole  mobility  is  independent  of 
the  doping  concentration  at  that  position.  If  you  set  MODVNP  to  1,  the 
dependence  is  given  by  the  other  parameters.  (This  option  has  not 
proven  useful  for  several  reasons  and  will  not  be  detailed  here.) 

Lines  8  and  9  include  the  parameters  determining  the  electron  mobil¬ 
ity,  The  equations  are  in  the  same  form  as  equations  (11)  to  (13): 

=  MUM  (1  -  C5-E)  ,  E  <  C7  ,  (14) 

=  C6-E-^/2(i  _  ciO-E-3/2)  ^  E  >  C7  ,  (15) 

//p  =  C16  +  C17-E  +  C18-E2  +  C19-E3  ,  C7  <  E  <  C20  .  (16) 

Line  10  has  been  covered  above. 


Line  11.  You  would  use  line  11  primarily  for  secondary  processes  in 
gases.  However,  for  semiconductors,  JMO  and  JPO  are  the  boundary 
values  for  the  electron  current  density  at  the  negative  electrode  and 
the  hole  current  density  at  the  positive  electrode,  respectively.  These 
are  used  primarily  under  reverse  bias. 

Line  12.  You  enter  voltage  parameters  on  line  12.  The  external  circuit 
used  in  DIODE  is  shown  in  figure  1  (many  of  the  labels  on  figure  1 
correspond  to  parameters  on  line  12;  see  table  1).  USTAT  is  the  initial 
source  voltage  and  V  is  the  initial  voltage  across  the  diode.  At  the  time 
Tl,  a  voltage  increment  DU  is  added  to  USTAT.  You  may  also  add  a 
sinusoidal  voltage  to  USTAT,  so  that  the  total  applied  voltage  is 

USTAT  +  VO  sin  OMEGA  t  ,  (17) 

where  t  is  the  problem  time.  You  may  use  this  option  to  apply  a 
constant  dV/dt,  since  for  OMEGA  t «  1,  dV/dt  is  essentially  equal  to 
VO  OMEGA. 


Figure  1.  Circuit  model 
used  in  DIODE. 
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Line  13.  You  may  apply  other  frequencies,  usually  harmonics,  by  using 
line  13,  PHIl  and  PHI2  are  the  phase  lags,  with  reference  to  OMEGA, 
for  OMEG  A1  and  OMEG  A2,  respectively.  IFOUR  may  be  set  at  one  to 
calculate  the  efficiency  of  oscillations  at  the  frequency  OMEGA  and  its 
second  harmonic.  If  IFOUR  =  1,  the  computer  calculates  a  cosine 
instead  of  sine  in  equation  (17).  IFOUR  is  in  column  80  of  line  13. 

Line  14.  Enter  geometrical  information  for  the  diode  on  line  14.  P  is  the 
gas  pressure  in  torr  and  should  be  set  equal  to  1.0  for  semiconductors. 
D  is  the  device  width  in  centimeters.  S  is  the  device  area  in  square 
centimeters.  T  allows  a  time  other  than  zero  to  be  the  initial  time  (this 
is  useful  in  extending  a  previous  nm  or  changing  the  phase  of  a 
sinusoidal  voltage).  DVDT  is  the  initial  dV/dt  used  in  the  start  of  the 
problem;  it  determines  the  displacement  current  at  zero  time,  but 
setting  dV/dt = 0  or  1  x  10^^,  for  example,  makes  little  difference  in  most 
problems.  JMAX  is  used  to  stop  calculations  when  unreasonably  high 
current  densities  are  calculated,  as  occurs  for  instabilities.  When  JA^ 
(see  sect.  4)  exceeds  JMAX,  the  calculations  are  halted  and  the  results 
of  that  time  step  are  printed  out. 

Lines  15  to  17.  Lines  15  to  17  are  used  to  describe  the  external  circuit. 
The  format  for  each  line  is  identical,  and  each  R,  L,  and  C  is  shown  in 
figure  1.  VC  and  DVC  are  the  initial  voltage  and  dV/dt  on  each  ca¬ 
pacitor,  used  mainly  for  extensions  of  previous  runs.  C  in  figure  1  is 
given  by  C  on  line  17.  Values  of  C  on  lines  15  and  16  are  ignored.  (The 
RLC  circuit  2  has  proven  useful  to  simulate  a  resistive  load,  but  not  as 
a  resonant  circuit,  because  transients  last  too  long  for  the  calculation 
to  be  co.'st  effective.  We  plan  to  modify  the  external  circuit  portion  of 
the  program  by  removing  RLC  circuit  3,  removing  L  from  the  supply 
circuit,  and  adding  LI  and  R1  to  the  diode  leg;  see  fig.  2.)  To  remove 
C2  and  C3  from  their  respective  legs  (see  fig.  1),  set  them  greater  than 
lx  10^0  F. 

Line  1 8.  Miscellaneous  decimal  data  are  entered  on  line  18.  STEPFA,  or 
F,  indicates  the  fraction  of  one  distance  step.  Ax,  that  the  electron  (or 
hole)  with  the  highest  velocity,  travels  in  one  time  step.  At.  This 
follows  from  the  calculation  of  Af  from 

At  =  F(Ax)/v^„^  .  (18) 

The  factor  f  may  be  used  to  ensure  that  the  calculation  time  step  is  less 
than  the  dielectric  relaxation  time  [7].  To  ensure  symmetric  diffusion, 
you  should  set  f  =  1/2  [8]. 

DJMAXM  is  a  check  on  too  rapid  a  current  increase  in  one  time  step. 
If  the  new  current  exceeds  the  old  current  by  a  factor  greater  than 
DJMAXM,  DT  is  halved  and  the  computation  repeated.  The  problem 
is  terminated  after  three  failures.  DIELK  is  the  dielectric  constant  of 
the  semiconductor.  DENSPL  is  the  carrier  density  of  holes,  in  cm"^  at 


Figure  2.  Proposed 
circuit  model  for 
DIODE,  adding 
capability  for  includ* 
ing  series  resistance 
and  inductance. 


the  right  boundary,  and  DENSMI  is  that  of  electrons  at  the  left 
boundary.  In  practice,  these  are  used  for  forward-biased  diodes  only. 
DENS2  is  the  square  of  the  intrinsic  density  at  300  K.  REC  is  the 
recombination  rate  for  the  semiconductor.  The  lifetime  in  intrinsic 
material  is  given  by  (2  REC)"^.  The  minority  lifetimes  of  electrons,  t„, 
and  holes,  tp,  are  given  by 

t„  =  n,/(REC)  Po  ;  tp  =  m,/(REC)  ,  (19) 

where  m,.  is  the  intrinsic  density  and  and  are  the  majority  carrier 
densities.  Finally,  EG  AP  is  the  energy  band  gap  in  the  semiconductor 
(see  eq  (8)). 

Line  19.  Line  19  determines  the  parameters  to  be  printed  out  as  a 
function  of  distance,  x,  at  selected  time  intervals.  The  first  two  param¬ 
eters  are  permanently  selected  to  be  x  and  E,  the  electric  field.  Users 
can  choose  six  additional  parameters  to  be  printed  out  from  the  list  of 
20  given  in  table  2. 

Line  20.  Line  20  includes  8  integers.  M  is  the  number  of  space  intervals. 
Ax  =  d/M.  M  +  1  values  of  each  distribution  parameter  must  be  in¬ 
cluded  in  the  input  data  (see  below.  Other  lines). 

PRTFRQ  is  the  print  frequency:  When  PRTFRQ  is  set  to  1,  M + 1  values 
are  printed  out;  for  PRTFRQ  =  2,  M/2  + 1  values  are  printed,  and  so 
on  for  higher  values  of  PRTFRQ. 

N3S  is  the  number  of  triples  to  be  included  on  line  21  (see  discussion 
of  line  21). 

MODP  is  the  mode  of  printout:  If  MODP  is  set  to  1,  the  printout  is 
according  to  time;  if  MODP  =  2,  the  printout  is  according  to  current 
(see  discussion  of  line  21). 

MODFCH  determines  whether  fixed  charges,  as  for  doped  semicon¬ 
ductors,  are  to  be  included.  If  MODFCH  =  0  or  1,  no  fixed  charges  are 
included;  if  MODFCH  >  1,  you  must  supply  fixed  charges  to  simulate 
the  device  doping. 
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Table  2.  Parameters 
inTPRINTLIST 
(line  21)  that  cause 
optional  output  data 
^o  be  printed 

Note:  A  maximum 
of  six  may  be 
chosen. 


Option 

No. 

Name 

Meaning 

1 

ALPHA 

Electron  ionization  coefficient,  a 

2 

BETA 

Hole  ionization  coefficient,  P 

3 

UNIVD 

Excitation  coefficient,  5 

4 

VM 

Electron  velocity 

5 

VP 

Hole  velocity 

6 

ALPHA*!- 

«/-(X) 

7 

BETA*J+ 

jS/+(X) 

8 

NP-NM  +  DN 

N+(X)  -  N_(X)  +  doping  density 

9 

RECOMRATE 

Recombination  rate 

10 

JPDIF 

Hole  diffusion  current  density 

11 

JNDIF 

Electron  diffusion  current  density 

12 

J  DISPL 

Displacement  current  density 

13 

J  TOTAL 

Total  current  density 

14 

VSUM 

Voltage  distribution 

15 

DVXDT 

dV/dt 

16 

N+(X) 

Hole  number  density 

17 

N-(X) 

Electron  number  density 

18 

J+(X) 

Hole  current  density 

19 

J-(X) 

Electron  current  density 

20 

TEMP 

Temperature 

EXTRAP  allows  extrapolations  to  be  used  instead  of  the  boundary 
conditions  given  in  line  11  or  18.  If  EXTRAP  is  set  to  0,  the  conditions 
in  lines  1 1  and  18  are  used;  otherwise  the  extrapolations  are  used.  (The 
extrapolation  option  has  rarely  been  used.) 

MS  sets  a  limit  on  the  number  of  consecutive  time  steps  in  which  a 
negative  field  is  calculated.  Calculation  is  halted  when  MS  is  ex¬ 
ceeded,  and  an  error  message  is  given. 

Last,  IPOINT  gives  the  number  of  time  steps  to  be  used  in  the 
summary  printout.  See  section  4  on  output. 

Line  22.  On  line  21  you  give  instructions  for  the  intervals  between  the 
time  steps  for  which  you  want  data  printed  out.  These  instructions, 
called  the  TPRINT  LIST,  are  composed  of  "triples":  two  real  numbers 
(including  a  decimal  point)  and  one  integer  compose  one  triple.  (The 
number  of  triples,  N3S,  was  specified  in  line  20.) 

If  MODP  (also  in  line  20)  is  set  to  1,  the  results  of  the  run  are  printed 
out  according  to  time.  In  this  case,  each  triple  is  interpreted  as  follows: 
the  first  number  is  the  time  for  the  initial  printout  and  ends  on  space 
14  of  line  21.  The  second  number  is  the  time  interval  between  subse¬ 
quent  printouts  and  ends  on  space  28.  The  integer  gives  the  number  of 
time  intervals  for  that  triple  and  ends  on  space  31.  (The  second  triple, 
if  used,  gives  the  same  information,  but  the  numbers  end  in  spaces  50, 
64,  and  67.) 


\ 
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If  MODP  =  2,  the  results  are  printed  out  according  to  current  density. 
In  this  case,  the  first  number  in  the  triple  gives  the  current  density  for 
the  first  printout,  the  second  number  gives  the  multiplier  for  subse¬ 
quent  printouts,  and  the  integer  gives  the  total  number  of  printouts. 
(The  spacing  of  the  fields  is  the  same  as  for  MODP  =  1.)  If  MODP  =  2, 
then  JMAX  on  line  14  gives  the  maximum  problem  time  in  seconds. 
Additional  lines  may  be  used  if  more  than  two  triples  are  desired. 

Other  lines.  Following  TPRINT  LIST  on  line  21,  you  can  use  a  variable 
number  of  lines  to  simulate  the  desired  initial  distributions  of  holes, 
electrons,  doping  levels,  and,  if  MODT  =  2  (line  4),  temperature.  You 
may  use  M  + 1  numbers,  5  to  a  line,  to  the  various  densities,  given  in 
units  of  cm"^.  The  first  density  number  must  end  on  space  14,  and  other 
numbers  on  multiples  of  14.  Enter  the  hole  distribution  first;  the 
electron  distribution  follows,  starting  with  a  new  line.  The  net  doping 
concentration  follows;  a  negative  sign  indicates  acceptors  or  p-type 
material.  Enter  the  initial  temperatures  in  kelvins. 

Ending  line.  Following  the  distributions  is  an  ending  line.  If  you  set 
ENDING  to  <10,000,  the  program  executes  another  nm.  If  ENDING 
exceeds  10,000  but  is  less  than  32,767,  the  program  stops. 


4.  Output 


The  output  of  DIODE  consists  of  a  printout  with  three  parts:  the  input 
data,  the  distribution  data  at  selected  times,  and  the  temporal  results. 
Appendix  B  is  a  portion  of  the  program  output  for  a  sample  problem. 

The  initial  page  of  the  input  data  is  headed  by  the  title  line,  reproduced 
from  the  input.  Each  of  the  input  lines  is  then  listed  in  a  format  similar 
to  the  input  and  identified  by  the  labels  shown  in  table  1.  (There  is  one 
exception  in  the  order  of  lines:  line  20  of  table  1  precedes  line  19.)  The 
TPRINT  LIST  (line  21)  prints  out  the  times  as  determined  by  the 
triples.  Following  the  TPRINT  LIST,  the  values  of  A6,  B6,  D6,  and  C9 
to  C19  are  printed  out.  (These  constants  are  explained  following  eq 
(10)  and  also  in  eq  (13).) 

The  second  page  of  the  printout  of  the  input  data  is  headed  by  the 
initial  time  and  voltage.  Then  the  distributions  follow  in  columns 
headed  by  M,  N+,  N-,  DN,  and  TEMP  K.  The  grid  points  range  from 
0  to  M,  for  M  +  1  values.  N+,  N-,  DN,  and  TEMP  K  are  the  input 
densities  of  holes,  electrons,  and  fixed  charges,  and  the  temperature, 
respectively. 

The  distribution  data  then  follow — one  page  for  each  time  listed  on  the 
TPRINT  LIST.  There  are  four  heading  lines  at  the  top  of  the  page.  Line 
1  includes  T  (time),  V  (diode  voltage),  U  (applied  voltage),  J+AVG 
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(average  hole  current  density),  J-AVG  (average  electron  current 
density),  J  AVG  (J+AVG  +  J-AVG),  and  GAMMA  (used  only  for  gas 
excitations). 

Line  2  includes  DT,  the  time  step;  V+MAX,  the  maximum  hole  velocity 
for  that  time  step;  V-M AX,  the  maximum  electron  velocity;  I  EXT,  the 
maximum  current  through  resistor  R;  I INT,  the  maximum  current 
through  the  diode;  J  INT,  the  average  current  density  (including  the 
displacement  and  diffusion  current  densities);  and  DV/DT,  the  rate  of 
voltage  rise  across  the  diode. 

Line  3  includes  D*ALPHAAVG,  ttie  integral  of  a  dx;  A*J-AVG,  the 
average  value  of  a/_;  D*BETAAVG,  the  integral  oihdx;  B^J+AVG,  the 
average  value  of  bj+;  RECOMRATAVG,  the  average  rate  of  electron- 
hole  recombination;  123,  the  sum  of  currents  through  circuits  2  and  3; 
and  DI/DT,  the  rate  of  current  rise  through  the  diode. 

Line  4  includes  J+DIF,  the  hole  diffusion  current  density  (calculated 
from  the  Einstein  equation  and  dp/dx  and  averaged);  J-DIF,  the  elec¬ 
tron  diffusion  current  density;  J  DIF,  the  total  diffusion  current  den¬ 
sity;  J  DISPL,  the  displacement  current  (calculated  from  dEfdT);  AVE 
TEMP  K,  the  average  temperature  across  the  diode;  TEM  RATE,  the 
rate  of  temperature  increase,  in  kelvins  per  nanosecond;  and  DENS2T, 
the  square  of  the  intrinsic  density  at  the  average  temperature. 

The  main  part  of  the  distribution  data  at  specified  times  is  arranged  in 
columns,  where  the  number  of  rows = (M  +  1)/PRTFRQ .  The  first  two 
columns  are  fixed  and  are  headed  X  (distance)  and  E(X)  (field).  The 
distance  is  given  in  centimeters  and  the  field  in  volts  per  centimeter. 
The  other  six  columns  are  chosen  by  the  user  from  the  options  shown 
in  table  2. 

The  page  following  the  distribution  data  is  headed  by  "capacitor 
voltages  at  last  successful  printout."  It  is  followed  by  the  voltages  and 
dV/df  for  the  capacitors  used  in  circuits  2  and  3.  The  final  page  (or 
pages)  of  the  printout  consists  of  six  columns  of  data  printed  accord¬ 
ing  to  time.  The  number  of  times  printed  is  the  value  of  IPOINT  (line 
20);  the  times  are  equally  spaced  between  the  initial  time  and  the 
maximum  time,  as  selected  by  the  TPRINT  LIST.  The  first  column  is 
headed  by  T  and  lists  the  time.  The  second  column  is  headed  by  V  and 
lists  the  voltages  for  each  time.  The  third  column  is  headed  by  I  and 
lists  the  diode  total  current.  The  fourth  is  headed  by  JAVG  and  is  the 
sum  of  the  electron  and  hole  current  densities.  The  fifth  column  is 
headed  by  II  and  lists  the  current  through  the  external  circuit  (see  fig. 
1).  The  last  column  is  headed  by  12  and  lists  the  current  through  circuit 
2.  This  completes  the  printout  for  the  program  DIODE. 


5. 


Material  Parameters 


The  first  10  lines  and  line  18  of  the  input  data  (see  table  1)  contain  the 
material  parameters  for  the  semiconductor  of  interest.  Table  3  con¬ 
tains  the  input  avalanche  parameters  for  silicon  and  gallium  arsenide 
which  have  been  most  commonly  used.  Table  4  contains  the  param¬ 
eters  used  for  the  variation  of  mobility  with  field.  For  silicon,  param¬ 
eters  are  given  for  high-mobility  (ideal)  and  low-mobility  (practical) 
material.  Other  material  parameters  (lines  4, 5,  and  18)  are  given  in 
table  5.  The  temperature  parameters  are  unknown  for  GaAs,  and  so 
the  silicon  values  are  used. 

From  equation  (11),  the  hole  velocity,  Vp,  is  given  by 

Vp(E)  =  MUP  •  £[1  -  (C1)E]  ,  E  <  C3  .  (20) 

Thus  the  hole  velocity  is  parabolic  with  a  maximum  velocity,  = 
MUP/4  •  Cl,  at  the  field  E^  =  1  /2(C1).  These  two  equations  for  Cl  are 
usually  not  compatible,  and  therefore  C3  must  be  chosen  less  thai  i  E^^^ 
and  equation  (13)  used  to  fit  the  experimental  curves.  Generally,  Vp 
calculated  from  equation  (12)  will  exceed  C4  where  E  is  slightly 
more  than  C7.  Tlie  same  considerations  apply  to  choosing  constants  in 
equations  (14)  to  (16)  to  fit  the  experimental  electron  velocity  versus 
field  curves. 

For  fitting  the  electron  velocity  versus  field  curves  for  GaAs,  or  other 
negative  differential  velocity  materials,  the  constant  CIO  in  equation 
(12)  must  be  chosen  to  be  negative.  Then  equation  (12)  has  a  minimum 
velocity,  at  the  field  E(v^jJ.  Then  the  fitting  equations  are 

EiVmin)  =  ~(2  •  C10)2/3  ;  v„,i„  =  (3/2)C6[E(i;„,-„)]i/2  .  (21) 

Again,  trial  and  error  will  be  required  to  obtain  the  best  fit  to  experi¬ 
mental  data. 

The  dependence  of  the  low  field  mobilities  of  the  electrons  and  holes 
upon  temperature  is  given  in  table  4  as  inverse  to  the  2.5  and  2.7  power, 
respectively.  However,  the  saturation  velocities  vary  approximately 
inversely  with  temperature.  At  present,  the  saturation  velocities  (C4 
and  C8  in  table  4)  are  increased  at  higher  ‘^emperatures  to  give  rough 
agreement  to  this  variation,  since  the  saturation  velocity  is  otherwise 
too  low.  A  planned  modification  to  the  program  will  give  a  single 
expression  [9]  that  gives  the  mobility  variation  with  field,  doping,  and 
temperature. 
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Table  3.  Parameter 
used  for  avalanche 
calculations 


Table  4.  Parameters 
used  for  mobility 
calculations 


Variable 

Value  for 
silicon 

Value  for 
GaAs 

MODA 

1 

1 

A1 

7.03  X 105 

5.05x106 

A2 

1.231  X 106 

1.95x106 

A3 

7.03x105 

1.52x106 

A4 

1.231  X 106 

1.60x106 

A5 

1.0  X  106 

2.86x105 

A6* 

1.231  X 106 

1.607x106 

MODS 

1 

1 

B1 

1.582x106 

2.05x107 

B2 

2.036x106 

2.49  X 106 

B3 

6.71  X 105 

1.49  X 106 

B4 

1.693x106 

1.75  X 106 

B5 

4.0x105 

2.82  X 105 

B6* 

1.693x106 

1.75x106 

CAT 

2x10-5 

2x10-3 

CBT 

2x10-3 

2  X 10-3 

*Computed 

Variable 

Value  for 
low-mobility 
silicon 

Value  for 
high-mobility 
silicon 

Value  for 
GaAs 

MODVP 

1 

1 

0 

MUP 

4.0x102 

6.0x102 

4.0x102 

Cl 

2.0x10-5 

4.0x10-5 

1.0x10-5 

C2 

1.6x104 

1.7x104 

4.47  xl04 

C3 

2.0x104 

9.0  X 103 

4.0x104 

C4 

1.0xl07 

1.0  X 107 

1.0xl07 

MODVM 

1 

1 

1 

MUM 

1.2x103 

1.5x103 

8.5x103 

C5 

4.4x10-5 

4.9x10-5 

1.0  X  1(H 

C6 

4.4x104 

4.4x104 

1.0xl04 

C7 

7.0x103 

1.0xl04 

5.5x103 

C8 

1.0xl07 

1.0  X 107 

1.0  X 107 

C9 

-1.55  xl07 

-1.55  X 107 

-4.77  xl06 

CIO 

-1.30  xl06 

-1.3x106 

2.4x107 

C15 

1.0xl05 

1.2  X 105 

5.0x104 

C20 

2.0x104 

2.0x104 

5.0x104 

Computed  parameters 

Cll 

2.90x106 

1.78  X 106 

0 

C12 

1.11  X 102 

2.05x102 

0 

C13 

-8.35x10-4 

-2.16x10-3 

0 

C14 

1.88  xir^ 

7.37x10-^ 

0 

C16 

1.97  xl06 

1.34  xl07 

2.63  xl07 

C17 

6.28  xl02 

-1.42  xl03 

-1.08  X 103 

C18 

-1.01  X 10-2 

1.09x10-4 

2.27x10-2 

C19 

-1.76x10-7 

-2.44  xl06 

-1.68x10-7 

Temperature  parameters 

PWRP 

2.7 

2.7 

2.1 

PWRN 

2.5 

2.5 

1.0 
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Table  5.  other  _  ...  „  Material 

material  parameters  Units  ^ 


DENSITY 

Density 

g/cm^ 

2.328 

5.32 

SPECHT 

Specific  heat 

J/g'C 

0.70 

0.35 

THRMCND 

Tbermal  conductivity 

W/cm®C 

1.45 

0.46 

DIELK 

Relative  dielectric  constant 

— 

13.1 

DENS2 

Intrinsic  density  squared 

cm"^ 

2x10^° 

3.2x10 

EGAP 

Energy  gap 

eV 

1.43 

CAT 

Proportional  change  in 
temperature  of  a 

2.0x10"* 

2.0x10 

•X 

CBT 

Proportional  change  in 
temperature  of  ^ 

2.0x10 

2.0x10 

PWRM 

Mobility  dependence  on 
power  of  temperature 
for  electrotrs 

2.5 

2.5 

PWRP 

Mobility  dependence  on 
power  of  temperature 
for  holes 

2.7 

2.7 

6.  Using  the  Interactive  Version  of  DIODE  on  an  IBM 
Terminal 

This  section  describes  how  to  use  a  dedicated  IBM  terminal  connected 
directly  to  the  HDL  IBM  3090  mainframe  to  access  the  DIODE  pro¬ 
gram.  Implementations  of  DIODE  on  other  hosts  will  differ  in  detail 
from  those  given  here  (particularly  instructions  for  logging  on,  etc). 
However,  certain  aspects  of  this  description  may  apply  to  any  imple¬ 
mentation  of  the  program. 

This  section  covers  procedures  for  logging  on,  defining  the  input  data, 
and  getting  the  output.  It  should  be  noted  that  this  document  gives 
only  a  cursory  description  of  the  IBM  3090's  menu-driven  work 
environment  (ISPF — interactive  system  productivity  facility)  and  of 
its  editor,  which  are  necessary  to  execute  the  program  and  specify  the 
input.  A  fuller  description  of  each  can  be  foimd  in  standard  IBM 
manuals  available  through  the  information  management  group  at 
HDL  or  by  contacting  Steven  Kaplan  (ext.  41403). 

6.1  Getting  Started 

These  instructions  assume  a  user  with  an  IBM  terminal  that  is  directly 
wired  to  the  HDL  IBM  3090.  The  first  step  in  using  the  DIODE  code  is 
to  log  onto  the  IBM  time  sharing  option  (TSO)  and  then  get  into  the 
menu-driven  mode  called  ISPF.  The  interactive  version  of  DIODE 
currently  resides  only  in  the  account  HK1005.  The  instructions  given 
here  are  easily  generalized  to  include  the  code's  use  on  other  accounts. 
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The  following  steps  are  required  to  access  the  main  ISPF  menu: 

(1)  Turn  on  the  terminal.  Make  sure  the  key  (if  there  is  one)  is  in  the 
horizontal  position.  Hit  RESET  if  the  keyboard  input  is  not  being 
accepted.  (If  there  is  still  a  problem,  contact  user  assistance,  ext.  42940.) 

(2)  Type  "LCX30N"  and  then  press  ENTER. 

(3)  When  prompted,  t)q)e  the  userid  (HK1005)  and  then  press  ENTER. 

(4)  When  prompted,  type  the  password  (see  Kaplan)  and  then  press 
ENTER. 

You  should  now  be  in  TSO.  After  some  system  information  and 
messages  appear  on  the  screen,  the  prompt  "READY"  should  appear 
at  the  left  end  of  the  current  line. 

(5)  At  the  READY  prompt,  type  "SPF". 

You  should  now  see  the  main  ISPF  menu,  which  reads  as  follows: 

ISPF  PRIMARY  MENU 


OPTION  ===> 

0  ISPFPARMS 

1  BROWSE 

2  EDIT 

3  UTILITIES 

4  FOREGROUND 

5  BATCH 

6  COMMAND 

7  dialoc;test 

8  LM  UTILITIES 
C  CHANGES 

S  SYSTEM  APPLICATIONS 
T  TUTORIAL 
U  USER  APPLICATIONS 
X  EXIT 

The  SPF  menu  has  14  options.  At  the  arrow  prompt  on  the  top  line, 
enter  tiie  number  or  letter  corresponding  to  the  desired  option.  The 
three  options  that  are  particularly  important  for  using  the  DIODE 
program  are  listed  below  with  brief  descriptions  of  their  respective 
functions. 


Option 

Mode 

Description 

1 

Browse 

Allows  user  to  view  a  dataset 

2 

Edit 

Allows  user  to  view  and  change  a  dataset 

6 

Command 

Allows  user  to  execute  a  TSO  command 

The  function  (PF)  keys  at  the  top  of  your  keyboard  perform  operations 
such  as  scrolling  through  datasets  and  moving  between  SPF  modes. 
The  PF3  key  allows  one  to  exit  from  the  different  SPF  modes  to  the 
main  menu.  The  uses  of  the  most  important  function  keys  are  as 
follows: 

PFl  Help  fimction 
PF3  Exit  current  mode 
PF7  Scroll  up  (edit/browse) 

PF8  Scroll  down  (edit/browse) 

PFIO  Scroll  left  (edit/browse) 

PFll  Scroll  right  (edit/browse) 

The  remainder  of  this  section  describes  using  the  SPF  options  to  run 
the  DIODE  program,  and  view  or  change  its  input/output. 

6.2  How  to  Edit/Browse  Datasets 

To  create  or  change  an  input  file  or  examine  an  output  file,  you  need 
to  use  tht  edit  or  the  browse  (no  changes  allowed)  mode.  The  follow¬ 
ing  describes  how  to  access  the  edit  mode.  (To  use  the  browse  option, 
follow  the  same  procedure,  except  enter  option  1  at  the  top  line  of  the 
SPF  main  menu.)  To  edit: 

Enter  option  2  at  the  arrow  prompt  on  the  top  line  of  the  SPF  main 
menu.  You  must  next  specify  the  file  to  be  edited.  A  dataset  panel 
appears  as  follows: 

ISPF  LIBRARY 

PROJECT  -=>  {enter  userid} 

GROUP  ==>  {enter  filename) 

TYPE  ==>  {enter  file  type} 

The  datasets  that  are  relevant  to  DIODE  are  the  following: 

INDIOD  DATA  :  input  dataset 

OUTDIOD  DATA  ;  output  dataset 

SKDIOD.FORT  :  source  code  dataset 

DIOD.CLIST  :  command  list  to  nm  the  program 
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So,  for  example,  to  edit  the  input  dataset,  the  panel  should  look  as 
follows: 


ISPF  LIBRARY 

PROJECT  ==> 

HK1005 

GROUP  ==> 

INDIOD 

TYPE  ==> 

DATA 

After  the  panel  is  set  up  to  indicate  the  appropriate  dataset,  hit 
ENTER.  The  dataset  contents  should  appear  on  the  screen.  You  can 
scroll  through  the  dataset,  and  change  its  contents  (if  in  edit  mode). 
When  editing,  be  sure  to  save  your  changes  by  typing  "SAVE"  at  the 
command  line  prompt  at  the  top  of  the  screen.  To  exit,  hit  the  PF3  key 
(this  should  also  save  the  dataset,  but  use  the  save  command  to  be 
sure). 

Note:  A  detailed  description  of  the  editor  commands,  the  function 
keys,  or  any  other  TSO/SPF  functions  can  be  obtained  by  contacting 
Steven  Kaplan  at  extension  41403. 

.3  How  to  Run  the  Program 

In  the  IBM  account  HK1005  there  is  a  command  list  (clist)  member 
called  DIOD(RD)  which  automatically  allocates  the  necessary  datasets, 
executes  the  program,  and  prints  the  output.  Variations  in  use  of 
DIODE,  such  as  incorporating  different  files  for  input  and  output,  are 
easily  accommodated  by  contacting  Steve  Kaplan  (ext.  41403).  To  run 
DIODE  using  the  existing  clist,  follow  these  steps: 

(1)  Starting  from  the  SPF  main  menu  (see  sect.  6.1),  get  into  the  TSO 
command  mode  by  entering  option  6  at  the  arrow  prompt  on  the  top 
line. 

(2)  In  TSO  an  arrow  prompt  will  also  appear.  At  this  prompt,  enter  the 
following: 

EXEC  DIOD(RD) 

where  EXEC  is  the  command  to  execute  a  clist  and  DIOD(RD)  is  the 
clist  to  be  executed. 

(3)  When  the  run  is  completed,  three  asterisks  {***)  will  appear.  Hit  the 
PF3  key  twice  to  return  to  the  SPF  main  menu. 


6.4.  Ending  the  Session 

Starting  from  the  main  SPF  menu: 

(1)  Type  "x"  or  "end"  at  the  arrow  prompt,  then  hit  ENTER. 

The  cursor  should  be  at  a  field  marked  "PROCESS  OPTION"  in  the 
subsequent  panel.  Enter  one  of  the  options  listed  in  this  panel;  in 
general,  either  "D"  (to  delete  the  dataset  containing  the  session  input) 
or  "K"  (to  keep  the  session  input — this  is  usually  unnecessary). 

(2)  Y ou  should  now  see  the  TSO  "READY"  prompt.  Log  off  the  system  by 
typing  "LO"  and  then  ENTER. 

7.  The  Program 

The  Fortran  77  program  is  divided  into  the  main  program  with  263 
source  statements  and  12  subroutines.  The  subroutines  are  listed  in 
table  6  with  the  purpose  and  number  of  statements  for  each.  Some 
statements  consist  of  multiple  lines. 

The  complete  listing  of  the  program  is  given  in  appendix  A.  A  portion 
of  the  printout  for  a  sample  problem  is  given  in  appendix  B. 

If  non-HDL  readers  are  interested  in  implementing  DIODE  on  other 
computers,  they  should  contact  the  authors  at  HDL  for  a  tape. 


Table  6.  Subroutines 
of  DIODE 


Title 

No.  of  source 
statements 

Purpose 

INPUT 

148 

Reads  and  writes  input  data 

INCON 

103 

Initial  calculations 

EXTCIR 

37 

Calculates  effect  of  external  circuit 

TEST 

72 

Tests  for  various  conditions 

OUTPUT 

87 

Calculates  and  prints  output  data 

ABD 

8 

Calculates  ionization  coefficients  (eq  (8),  (9)) 

UNIV 

24 

Further  ionization  calculations 

MUF 

11 

Calculates  mobilities  (eq  (10)-(15)) 

TESTSQ 

29 

Gas  excitation 

SIMQ 

53 

Further  gas  excitation 

EFFCY 

98 

Calculation  of  rf  efficiency 

BLKDT 

20 

Reserves  memory  blocks 

Statements  in  main  program; 

263 

Statements  in  subroutines: 

690 

Total 

953 

23 
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Appendix  A 


c 

c 


c 


2 


10 


C 

C 

C 


DIODE  SPACE  CHARGE  PROGRAM  FOR  GASES  AND  SEMICONDUCTORS  WITH  DRIFT,  A 


IONIZATION  AND  RECOMBINATION  AS  A  FUNCTION  OF  TEMPERATURE  A 

REAL  ICUR, ICURO, lEXT, I INT, I  INTO, I23,JAVG, JC, JCSAV,JDIF, JDISPL,JDIS  A 
ISP , J I  NT , JM , JMAX , JMO I F , JMSAV , JMO , JP , JPD I F , JPSAV , JPO ,MJDSP ,MJM,ltJMDF  A 
2,MJP,MJPDF,MTEM,MUM,MUMO,MUP,MUPO,NM,NMSAV,NP,NPP,NPSAV,NOM,NOP,L,  A 
3MUF,MUMF,MUPF  A 

DIMENSION  DVXDT(101),DNPDX(101),DNMOX(101),DIFP(101),DIFM(101),  A 
1GM(101),GP(101),JDISSP(101),OTEM(101)  A 

INTEGER  PRTFRQ, CHECK, ENDING, CHECK2,EXTRAP,CIRCT( 10), GHECK1  A 

COMMON  /ABC/  MODA,MODB,MODD,MODVP ,MODVM,MODVNP ,M0DVNM,M0DT,A1 ,A2,  A 
1A3,A4,A5,A6,B1,B2,B3,B4,B5,B6,01,D2,D3,D4,D5,06,C1,C2,C3,C4,C5,C6,  A 
2C7,C8,C9.C10,C11,C12,C13,C14,C15,C16,C17,C18,C19,C20,MUP,MUM,MUP0,  A 
3MUM0,AVP,AVM, NOP, N0M,VP1M,VM1M,CAT,CBT,PWRP.PWRM, TEMP  A 

COMMON  /CIRT/  USTAT,U,V,DU,V0,V1 ,V2,0MEGA,0MEGA1 ,0MEGA2,PH 1 1 ,PHI2,  A 
1RD,LD,R(3),L(3),C{3),CP,T1,IEFCY,VNEW  A 

COMMON  /DSTR/  NP (201 ) ,NM(201 ) ,DN(201 ) ,TEM(201 ) ,E(201 ) ,DE(201 ) ,VSUM  A 
1(201),VP(201) ,VM(201) ,JP(201),JM(201),JPDIF(201) ,JMDIF(201) ,JDISPL  A 
2(201),RCMBR(201) ,ALPHA(201) ,BETA(201),TPRINT(201)  A 

COMMON  /EFCY/  AMP, AMP2, PHASE, PJ,PV,PHASE2,PJ2,PV2, DTP, GNEG,RFP,EFF  A 
1 ,AMPVDC,AMPJDC,0CPWR,AMPV,AMPV2  A 

COMMON  /FUND/  M.MPI ,N3S, PRTFRQ, NXPRNT(6) ,MODP,MODFCH, EXTRAP, MS, IJ,  A 
1 IPOINT,FM,STEPFA,DX,DT,T,D,S,P,EL,EPSO,DIELK,DNSTY,SPCHT,THCND,DEN  A 
2SPL,DENSMI ,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI , GAMMA, TAUSEC  A 

COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT,l INT,I INT0,IEXT,I23,DIDT,  A 
10V0T,DVDT0,VC(3),VC0(3),DVC(3),0VC0(3),ICUR(3),  ICUR0|3),MJP,MJM,  A 
2MJPDF,MJMDF,MJ0?.P,MTEM,VPMAX,VMMAX,CST(10,3),DMTEM,D0IDT,DIDT0  A 

COMMON  /PLO/  TIME(401),VOLT(401),CURR(401),JAVG(401),CURR1(401),  A 
1CURR2(401),CURR3(401),CLTT(401) ,TAVR{401)  A 

COMMON  /SAV/  NPSAV(101),NMSAV(101),ESAV(101),VXSAV(101),VPSAV(101)  A 
1 ,VMSAV( 101 ) , JPSAV( 101 ) , JMSAV( 101 ) ,VSAV,TSAV, JCSAV,GAMSAV,TMPSAV  A 

A 

MUPF(Z)=MUF(Z,P,C1,C2,C3,C4,C9,MUP)  A 

MUMF(Z)=MUF(Z,P,C5,C6,C7,C8,C10,MUM)  A 

UNIVA{EA)=UNIV(P,M0DA,EA,A1,A2,A3,A4,A5,A6)  A 

UNIVB(EA)=UNIV(P,M0DB,EA,B1,B2,B3,B4,B5,B6)  A 

UNIVD(EA)=UNIV(P,M0DD,EA,01,D2,D3,D4,D5,D6)  A 

B0LTZ=8.6E-5  A 

ENDING=10  A 

CHECK  1=0  A 

CALL  INPUT(ENDING,CHECK1)  A 

VNEW=V 

IF  (CHECKI.EQ.I)  GO  TO  275  A 
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20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 
31B 
31C 
31D 
31E 
31F 

32 

33 
33A 
33B 

33C 


SET  COUNTERS 


A  34 


INSTEP=0 

DT=0. 

I  I NT=0 

DIDT=0 

DDIDT=0 

ISET=1 

INDEX=1 

M I NUS=MS 


A  35 
A  36 


A  37 
A  38 
A  39 
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MSO=MS  A  40 

CHECK=0  A  41 

CHECK2*0  A  42 

C 

C  COMPUTE  CONSTANTS  FROM  INPUT  A  43 

C 

CALL  INCON(CIRCT, ENDING)  A  44 

IT=0  A  45 

DTP=(TPRINT( IJ)-TPRINT(1))/(FL0AT(IP0INT)-1.)  A  46 

IF  ( INSTEP. EQ.O)  GO  TO  125  A  47 

C  START  OF  COMPUTATION  LOOP  A  48 

110  INSTEPsI  A  49 

C  COMPUTE  NEW  TIME  A  50 

DT=STEPFA*DX/VMAX  A  51 

GO  TO  123  A  52 

120  DT=DT/2.  A  53 

123  T=T+DT  A  54 

C  COMPUTE  U  A  55 

U=USTAT  A  56 

IF  (T.LT.T1)  GO  TO  125  A  57 

IF  (lEFCY.EQ.O)  GO  TO  124  A  58 

U=U+DU+V0*C0S(0MEGA*T)+V1*C0S(0MEGA1*T+PH1 1 )+V2*C0S(0MEGA2*T+PH 12)  A  59 

GO  TO  125  A  60 

124  0=U+DU+V0*SIN(0MEGA*(T-T1))  A  61 

125  CONTINUE  A  62 

DO  126  l=1,HP1  A  62A 

GP(I)=EL*NP(I)  A  62B 

126  GM(I)=EL*NM{I)  A  63 

C  CALCULATE  RECOMBINATION  RATE  A  64 

0ENS2T=DENS2» ( MTEM/300 . ) **3*EXP ( EGAP/BOLTZ* ( 1 . /300 . - 1 . /MTEM ) )  A  65 

DO  128  1=1, MP1  A  66 

NPP=NP(I)»NM(I)  A  67 

IF  {0ENS2T.NE.0,)  NPP=(NPP-DENS2T)/SQRT(DENS2T)  A  68 

128  RCMBR( l)=REC*EL*NPP  A  69 

IF  ( INSTEP. EQ.O)  GO  TO  169  A  70 

C  CALCULATION  OF  EXTERNAL  CIRCUIT  -  V.JEXT  A  70A 

CALL  EXTCIR  (ISET.CIRCT)  A  70B 

C  HOLE  CHARGE  DENSITY  A  71 

IF  (EXTRAP.EQ.1)  GO  TO  130  A  72 

GP(MP1)=DENSPL*EL  A  73 

GO  TO  140  A  74 

130  GP(MP1)=GP(MP1)+(2.*JP(MP1)-3.*JP(M)+JP(M-1))«DT/DX  A  75 

140  SUM=(UNIVD(E(MP1))*JM(MP1)+UNIV0(E(1))«JM(1))/2.  A  76 

DO  150  l=1,M  A  77 

GP( I )=GP( I )+DT*(ALPHA( I )*JM( I )+BETA{ I )»JP( I )+( JP( l+1)-JP( I ) )/DX)  A  78 

150  SUM=SUM+UNIVD(E( l))*JM{ I)  A  79 

C  COMPUTE  EFFECT  OF  PHOTONS  HITTING  CATHODE  A  80 

RFAC=0.  A  81 

IF  (TAUSEC.NE.O.)  RFAC=EXP(-DT/TAOSEC)  A  82 

GAMMA=RFAC*GAMMA+GAMSEC*DX*(1.-RFAC)*SUM  A  83 

C  ELECTRON  CHARGE  DENSITY  A  84 

DO  160  1=2, MP1  A  85 
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160  GM(I)=DT»(ALPHA(I)*JM(I)+BETA(I)*JP(I)-(JM(I)-JM(I-1))/0X)+GM(I)  A  86 

DO  162  1=1, MP1  A  87 

GP(I)=GP(I)-DT*RCMBR(I)  A  88 

162  GM(I)=GM(I)-DT«RCMBR(I)  A  89 

C  TEST  FOR  NEGATIVE  NUMBER  DENSITIES  A  90 

DO  168  1=1, MP1  A  91 

NP(I)=GP(I)/EL  A  92 

NM( l)=GM( l)/EL  A  93 

IF  {NP(l).GE.O.)  GO  TO  165  A  94 

GP(l)=0.  A  96 

NP(l)=0.  A  97 

165  IF  (NM(I).GE.O.)  GO  TO  168  A  98 

GM(l)=0.  A  99 

NM(l)=0.  A  100 

168  CONTINUE  A  102 

C  ELECTRIC  FIELD  A  103 

169  ESTAR=0.0  A  104 

E(1)=0.  A  105 

E(2)=GM(1)+GM(2)-GP(1)-GP(2)-EL*(DN(1)+DN(2))  A  105A 

DO  170  1=3, M  A  106 

170  E(l)=E(l-1)+(GM(l-2)+2.<»GM(l-1)+2.»GM(l)+GM(l+1)-GP(l-2)-2.»GP(l-1  A  107 

1)-2.»GP(l)-GP(l+1))/3.-EL*(DN(l-1)+DN(l))  A  107A 

E(MP1)=E(M)+GM(M)+GM(MP1)-GP(M)-GP(MP1)-EL*(DN(M)+DN(MP1))  A  107B 

DO  180  1=2, MP1  A  107C 

180  ESTAR=ESTAR+E(I-1)+E(I)  A  108 

E(1)=VNEW/D-CST(6,1)<*ESTAR/(2.«*FM)  ’  A  109 

E(1)=E(1)-CST(6,1)*((GM(1)-GP(1)-GM(MP1)+GP(MP1)-EL^*(DN(1)-0N(MP1)  A  110 
1))/(6.*FM))  A  111 

C  ADJUST  E  A  112 

DO  190  1=2, MP1  A  113 

E(I)=CST(6,1)*E(I)+E(1)  A  114 

190  VSUM( I )=VSUM( l-1)+0.5*(E( I )+E( l-1))«0X  A  115 

C  COMPUTE  FUNCTIONS  OF  ELECTRIC  FIELD  A  116 

IF  (INSTEP  .EQ.O)  GO  TO  198  A  117 

DO  195  1=1, MP1  A  117A 

DE(I)=E(I)-ESAV(I)  A  118 

JDISPL( I )=DIELK*EPSO*DE( I )/DT  A  119 

195  DVXDT(I)=(VSUM(I)-VXSAV(I))/DT  A  120 

C  ELECTRON  VELOCITY-VM  A  121 

198  DO  230  1=1, MP1  A  122 

AOP=E( I )/P  A  123 

IF  (MODVM.EQ.O)  GO  TO  200  A  124 

IF  (A0P.LE.C7.0R.A0P,GE.C20)  GO  TO  200  A  125 

VMO=C16+C17*AOP+C18*AOP*AOP+C19*AOP*AOP»AOP  A  126 

GO  TO  201  A  127 

200  VM0=E(l)/P»MUMF(E(l))  A  128 

201  IF  (MODVNM.EQ.O)  GO  TO  205  A  129 

ADN=ABS(DN(I))  A  130 

DM0=(A0N/N0M)**AVM  A  131 

VM1=MUM0*A0P  A  132 

IF  (VM1.GT.VM1M)  VM1=VM1M  A  133 

VM(l)=(VM1+(VM0-VM1)/(1.+DM0))*(30O./TEM(l))»*PWRM  A  134 

GO  TO  210  A  135 
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205 

VH( 1 )=VM0*(300./TEM( 1 ) )»*PWRM 

A 

136 

C 

HOLE  VELOCITY-VP 

A 

139 

210 

IF  (MODVP.EQ.O)  GO  TO  220 

A 

140 

IF  (A0P.LE.C3.0R.A0P.GE.C15)  GO  TO  220 

A 

141 

VP0=C 1 1 +C1 2*A0P+C 1 3<*A0P*A0P+C 14«A0P*A0P*A0P 

A 

142 

GO  TO  221 

A 

143 

220 

VP0=MUPF(E(  1 )  )^*E{  1  )/P 

A 

144 

221 

IF  (MODVNP.EQ.O)  GO  TO  225 

A 

145 

DP0=(ADN/N0P)**AVP 

A 

146 

VP1=MUPO*AOP 

A 

147 

IF  (VPI.GT.VPIM)  VP1=VP1M 

A 

148 

VP( 1 )=(VP1+(VP0-VP1)/( 1 .+DP0) )»{300./TEM( I ) )«<*PWRP 

A 

149 

GO  TO  230 

A 

150 

225 

VP( 1 )=VP0*(300./TEM{ 1 ) )*»PWRP 

A 

151 

230 

CONTINUE 

A 

151A 

C 

COMPUTE  MAXIMUM  VELOCITIES 

A 

152 

VMMAX=VM(1) 

A 

152A 

VPMAX=VP(1) 

A 

152B 

DO  235  1=1, MP1 

A 

152C 

VMMAX=AMAX1 (VMMAX,VM( 1 ) ) 

A 

152D 

235 

VPMAX=AMAX1{VPMAX,VP{ 1 )) 

A 

152E 

VMAX=AMAX 1 ( VPMAX , VMMAX ) 

A 

152F 

C 

COMPUTE  DIFFUSION  CURRENT  DENSITIES 

A 

152G 

DO  255  1=1, MP1 

A 

152H 

DNMDX(  1  )  =  (NM(  l+1)-NM(  l-1))/(2.^»DX) 

A 

1521 

255 

DNPDX( l)=(NP( l+1)-NP( l-1))/(2.*DX) 

A 

152J 

DNMDX(1)=(NM(2)-NM(1))/DX 

A 

152K 

DNMD):(MP1)=(NM(MP l)-NM(M))/DX 

A 

152L 

DNPDX(1)=(NP(2)-NP(1))/0X 

A 

152M 

DNPDX(MP1)-(NP(MP1)-NP(M))/DX 

A 

152N 

DO  265  1=1, MP1 

A 

1520 

AOP=E( 1 )/P 

A 

152P 

DIFM( 1 )=0. 

A 

152Q 

IF  (AOP.NE.O.)  DIFM(I)=BOLTZ*TEM{I)*VM(I)/AOP 

A 

152R 

JMDIF( 1 )=-DIFM( 1 )*DNMOX( 1 )*EL 

A 

152S 

DIFP( 1 )=0, 

A 

152T 

IF  (AOP.NE.O.)  DIFP(l)=BOLTZ^^TEM(l)*VP(l)/AOP 

A 

152U 

265 

JPDIF( 1 )=DIFP( 1 )*DNPDX( 1 )*EL 

A 

152  V 

C 

CURRENT  OENSITIES-JP,JM 

A 

152W 

DO  240  1=1, MP1 

A 

153 

JP( 1 )=GP( 1 )*VP(  1 ) 

A 

154 

240 

JM(I)=GM(I)*VM(I) 

A 

155 

JP(MP1)=JP(MP1)+JP0 

A 

155A 

GP(MP1)=ABS(JP(MP1)/VP(MP1)) 

A 

155B 

NP(MP1)=GP(MP1)/EL 

A 

155C 

DENSIT=DENSMI 

A 

156 

IF  (EXTRAP. EQ.1)  DENSIT=NM(1) 

A 

156A 

C 

COMPUTE  INITIAL  GAMMA 

A 

157 

IF  ( INSTEP. EQ.1)  GO  TO  245 

A 

157  A 

IF  (GAMMA, GE.O.)  GO  TO  243 

A 

157B 

GAMMA=0 . 

A 

157C 

DO  242  1=2, M 

A 

1570 
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242  GAMMA=GAMMA+JM(I)^*UNIVD(E{I))  A  157E 

GAMMA= ( GAMMA+ ( UN  I VD ( E ( 1 ) ) * JM ( 1 ) +UN I VO ( E I MP 1 ) ) * JM ( MP 1 ) ) /2 . ) *GAMSEC*  A  1 57F 

1DX  A  157G 

243  IF  ( INSTEP. EQ.O)  GO  TO  250  A  157H 

245  JM{1)=JMO+GAMI*JP(1)+GAMMA+OENSlT*EL*VM(1)  A  158 

IF  (EXTRAP. EQ.O)  GO  TO  250  A  159 

JMn)-JM(1)+VM(1)*0T/0X*(2.*JM(1)-3.*JM(2)+JM(3))  A  159A 

250  GM(1)=ABS(JM(1)/VM(1))  A  159B 

NM(1)=GM(1)/EL  A  159C 

C  COMPUTE  AVERAGE  CURRENT  DENSITIES  A  160 

MJP=(JP(1)+JP(MP1))/2.  A  160A 

MJM=(JM(1)+JM(MP1))/2.  A  160B 

MJPDF=(JPDIF(1)+JPDIF(MP1))/2.  A  160C 

MJMDF=(JMDIF(1)+JM0IF(MP1))/2.  A  1600 

MJDSP=(JDISPL(1)+JDISPL(MP1))/2.  A  160E 

DO  400  1=2, M  A  161F 

MJP=MJP+JP( I )  A  160G 

MJM=MJM+JM( I )  A  160H 

MJPDF=MJPDF+JPD I F ( I )  A  160 1 

MJMDF=MJMDF+JMDIF( I )  A160J 

400  MJDSP=MJDSP+JDISPL( I )  A  160K 

MJP=MJP/FM  A  160L 

MJM=MJM/FM  A  160M 

MJPDF=MJPDF/FM  A  161 

MJMDF=MJM0F/FM  A  161A 

MJDSP=MJDSP/FM  A  161B 

JC=MJP+MJM  A  161C 

JOIF=MJMOF+MJPDF  A  1610 

JINT=JC+JDIF+MJDSP  A  161E 

I INT=ABS(S*JINT)  A161F 

IF(CIRCT(1).Eq.5)  IEXT=(U-V)/R(1)  A  161G 

I23=ICUR(2)+ICUR(3)  A  161H 

VNEW=V-RD»I INT-LD»DIDT 

IF  (DT.EQ.O)  DIDT=0.  A  1611 

DDIDT=0. 

IF(DT.NE.O.)  DIDT=( I INT-I INTO)/DT  A  161J 

I F ( DT . NE . 0 . )  DO  I DT= ( 0  I DT-D I DTO ) /DT 

IF  (MOOT. EQ.O.)  GO  TO  30  A  161K 

DO  410  1=1, MP1  A  162 

JDISSP( I )=JP( I )+JK; I )+JPDIF( I )+JMDIFI I )  A  163 

DTEM( I )=JDISSP( I )*E( I )*DT/(DNSTY*SPCHT)  A  164 

410  TEM(I)=TEM(I)+DTEM(I)  A  165 

MTEM=(TEM(1)+TEM(MP1))/2.  A  166 

DO  420  1=2, M  A  167 

420  MTEM=MTEM+TEM( I )  A  168 

MTEM=MTEM/FM  A  169 

DMTEM=0.  A  169A 

D2T=DT*1.E9 

IF  (DZT.NE.O.)  DMTEM=(MTEM-TMPSAV)/DZT  A  169B 

C  END  OF  COMPUTATION  LOOP  A  169C 

C  A52-91  INSERTED  BETWEEN  A169D  AND  170  A  169D 

C  TEST  FOR  STABILITY  AND  OTHER  CRITERIA  A  52 

30  CALL  TEST(MINUS, INSTEP, CHECK,CHECK2,ISET,MS0)  A  53 

GO  TO  (40,120,260),  ISET  A  54 


41 


Appendix  A 


40  CHECK2«0  A  55 

DO  45  l=1,MP1  A  56 

ALPHA(l)=UNIVA(E(l))*(1.-CAT*(TEM(l)-300.))  A  56A 

45  BETA(l)=UNIVB(E(l))*|1.-CBT*(TEM(l)-300.))  A  56B 

C  TEST  FOR  PRINTING  A  56C 

IF  (MODP.Eq.2)  GO  TO  60  A  57 

IF  (T.LT.(FLOAT( IT)*DTP+TPRINT(1)))  GO  TO  50  A  58 

IT=IT+1  A  59 

TIME(IT)=T  A  60 

V0LT(IT)=V  A  61 

CURR( IT)s| INT  A  62 

JAVG(IT)=JC  A  63 

comn(rT)=icuR(i)  a  64 

CURR2( IT)=ICUR(2)  A  65 

CURR3( IT)-ICUR(3)  A  66 

CLTT( IT)=DIOT 
TAVR( IT)»MTEM 

50  IF  (T.LT.TPRINT( INDEX))  GO  TO  100  A  67 

GO  TO  70  A  68 

60  IF  (JINT.LT.TPRINT( INDEX))  GO  TO  100  A  69 

G  CALCULATE  VALUES  TO  BE  PRINTED  AND  PRINT  THEM  A  70 

C  (72  TO  76  REMOVED)  A  71 

70  CALL  OUTPUT (DVXDT)  A  77 

INDEX=INDEX+1  A  78 

IF  ( INDEX. LE.IJ)  GO  TO  100  A  79 

WRITE  (6,360)  A  80 

WRITE  (6,350)  ( VC( I ) ,0VC( I ) , 1*1 ,3)  A  81 

IF( IEFCY.NE.1)  GO  TO  270  A  82 

CALL  EFFCY  A  83 

WRITE  (6,500)  AMP,AMPV,PJ,PV,PHASE,GNEG,RFP,AMPVDC,AMPJDC,  A  83A 

10CPWR,EFF  A  83B 

WRITE  (6,501)  AMP2,AMPV2,PJ2,PV2,PHASE2  A  83C 

GO  TO  270  A  84 

100  IF  (CHECK. EQ.O)  GO  TO  110  A  85 

CHECK*0  A  86 

MINUS=MINUS-1  A  87 

IF  (MINUS. GT.O)  GO  TO  110  A  88 

WRITE  (6,310)  A  89 

260  WRITE  (6,320)  A  90 

CALL  OUTPUT  (DVXDT)  A  91 

270  WRITE  (6,330)  A  170 

WRITE  (6,340)  (TIME( IK),VOLT( IK)-LD*CLTT( IK)-RD^»CURR( IK),CURR( IK) 

1,JAVG( IK),CURR1( IK),CURR2( IK),VOLT( IK]r,TAVR( IK) , IK=1, IPOINT)  A  172 

275  READ  (1,290)  ENDING  A  174 

IF  (ENDING. GE. 10000)  GO  TO  280  A  175 

WRITE  (6,300)  A  176 

GO  TO  10  A  177 

280  STOP  A  178 

C  A  179 

C  A  180 

290  FORMAT  (15)  A  181 

300  FORMAT  (1H1)  A  182 

310  FORMAT  (//51H  THERE  HAVE  BEEN  MS  CASES  OF  NEGATIVE  VELOCITIES. /I  A  183 

16H  END  OF  PROBLEM.)  A  184 
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320  FORMAT  (47H  PRINTOUT  AT  TIME  OF  ERROR  FOLLOWS  ON  NEXT  PAGE)  A  186 

330  FORMAT  ( 1H1 , 13X, 1HT, 14X, 1HV, 14X, 1H I , 1 1X,4HJAVG, 13X,2HI 1 , 13X,2H 12,  A  187 

18X,7HV(DLEG),7X,8HAVG.TEMP)  A  188 

340  FORMAT  (8(1PE15.6))  A  189 

350  FORMAT  (1P,6E20.8)  A  190 

360  FORMAT  ( 1H1 , 10X,46HCAPACIT0R  VOLTAGES  AT  LAST  SUCCESSFUL  PRINTOUT  A  191 

1//17X,3HVC1,16X,4HDVC1,17X,3HVC2,16X,4HDVC2,17X,3HVC3,16X,4HDVC3)  A  192 

C  FORMAT  (5E14.8)  A  193 

500  FORMAT  (//47H  THE  AMPLITUDE  OF  THE  FIRST  HARMONIC  OF  JAVG  IS,  A  194 

*1PE20.8/49H  THE  AMPLITUDE  OF  THE  FIRST  HARMONIC  OF  V(OLT)  IS,  A  195 

*1PE20.8/54H  THE  PHASE  OF  THE  FIRST  HARMONIC  OF  JAVG  IN  DEGREES  IS,  A  196 

*1PE20.8/51H  THE  PHASE  OF  THE  FIRST  HARMONIC  OF  V  IN  DEGREES  IS,  A  197 

*1PE20.8/43H  THE  PHASE  BETWEEN  V  AND  JAVG  IN  DEGREES  IS,  1PE20.8/  A  198 

*28H  THE  NEGATIVE  CONDUCTANCE  I S, 1PE20 .8/16H  THE  RF  POWER  IS,  A  199 

*1PE20.8/28H  THE  TIME  AVERAGE  VOLTAGE  IS, 1PE20.8/31H  THE  AVERAGE  CU  A  200 
*RRENT  DENSITY  IS,1PE20.8/15H  INPUT  POWER  IS,1PE20.8/18H  THE  EFFICI  A  201 
*ENCY  IS,1PE20.8)  A  202 

501  FORMAT  (//48H  THE  AMPLITUDE  OF  THE  SECOND  HARMONIC  OF  JAVG  IS,  A  203 

*1PE20.8/50H  THE  AMPLITUDE  OF  THE  SECOND  HARMONIC  OF  V(OLT)  IS,  A  204 

*1PE20.8/55H  THE  PHASE  OF  THE  SECOND  HARMONIC  OF  JAVG  IN  DEGREES  IS  A  205 
*,1PE20.8/52H  THE  PHASE  OF  THE  SECOND  HARMONIC  OF  V  IN  DEGREES  IS,  A  206 
*1PE20.8/43H  THE  PHASE  BETWEEN  V  AND  JAVG  IN  DEGREES  IS,1PE20.8)  A  207 

END  A  208- 

SUBROUTINE  INPUT(ENDING,CHECK1)  B  1 

INTEGER  PRTFRQ, EXTRAP, ENDING, CHECK1  B  2 

DIMENSION  TPRNTI(66),  TINC(66),  NINC(66)  B  3 

REAL  JPOIF,JMDIF,NP,NM,JP,JM,JDISPL,MJP',MJM,MJPDF,MJMDF,MJDSP  B  4 

REAL  MUM, MUP,JMAX,JPO,JC,MTEM,  L, I INT,JINT, I  INTO, ICUR, ICURO, lEXT,  B  5 
1 I23,MUP0,N0P,MUM0,N0M,JM0  B  6 

COMMON  /ABC/  M0DA,M0DB,M0DD,M0DVP,M0DVM,M0DVNP,M0DVNM,M0DT,A1 ,A2,  B  7 
1A3,A4,A5,A6,B1,B2,B3,B4,B5,B6,01,D2,D3,D4,D5,D6,C1,C2,C3,C4,C5,C6,  B  8 
2C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,C18,C19,C20,MUP,MUM,MUP0,  B  9 
3MUM0,AVP,AVM, NOP, N0M,VP1M,VM1M, CAT, CBT,PWRP,PWRM, TEMP  B  10 

COMMON  /CIRT/  USTAT,U,V,DU,V0,V1 ,V2,0MEGA,0MEGA1 ,0MEGA2,PH 1 1 ,PH 12,  B  11 
1RD,LD,R(3),L(3),C(3),CP,T1,IEFCY,VDD  B  12 

COMMON  /DSTR/  NP|201 ) ,NM(201 ) ,DN(201 ) ,TEM(201 ) ,E(201 ) ,DE(201 ) ,VSUM  B  13 
1(201),VP(201),VM(201),JP(201),JM(201),JPDIF(201),JMDIF(201),JDISPL  B  14 
2(201), RCMBR{201),ALPHA{201),BETA(201),TPRINT(201)  B  15 

COMMON  /FUND/  M,MP1 ,N3S, PRTFRQ, NXPRNT(6) ,MODP ,MODFCH,EXTRAP,MS, I J,  B  16 
1 IPOINT,FM,STEPFA,DX,DT,T,D,S,P, EL, EPSO, DIELK, DNSTY,SPCHT,THCND, DEN  B  17 
2SPL,DENSMI ,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI , GAMMA, TAUSEC  B  18 

COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT,IINT,IINTO,IEXT,I23,DIDT,  B  19 
1DVDT,DVDTO,VC(3),VCO(3),DVC(3),DVCO(3),ICUR(3),  ICUR0(3) ,MJP,MJM,  B  20 
2MJPDF,MJMDF,MJDSP,MTEM,VPMAX,VMMAX,CST(10,3),DMTEM,DDIDT,DIDT0  B  21 
READ  (1,120)  B  21A 

WRITE  (6,120)  B  22 

IF(ENDING.LE.3)  GO  TO  25  B  22A 

READ  (1,150)  M0DA,A1,A2,A3,A4,A5  B  23 

READ  (1,150)  M0DB,B1,B2,B3,B4,B5  B  24 

READ  (1,150)  M0DD,D1,D2,D3,D4,D5  B  25 

READ  (1,150)  MOOT, CAT, CBT,PWRM,PWRP, TEMP  B  25A 

READ  (1,130)  DNSTY,SPCHT,THCND  B  25B 

READ  (1,150)  M0DVP,MUP,C1,C2,C3,C4  B  26 

READ  (1,150)  M0DVNP,MUP0,AVP,N0P,VP1M  B  26A 


43 


Appendix  A 

READ  (1,150)  M0DVM,MUM,C5,C6,C7,C8  B  27 

READ  (1,150)  M0DVNM,MUM0,AVM,N0M,VM1M  B  27A 

READ  (1,130)  C9,C10,C15,C20  B  28 

READ  (1,130)  GAMSEC, TAUSEC, GAMMA, GAMI ,JM0,JP0  B  29 

READ  (1,130)  USTAT,DU, VO, OMEGA, T1,V  B  30 

READ  (1,500)  V1,0MEGA1,PHI1,V2,0MEGA2,PHI2,IEFCY  B  30A 

READ  (1,130)  P,D,S,T,DVDT,JMAX  B  31 

DO  10  IGIR=1,3  B  32 

10  READ  (1,130)  R(ICIR),L(1GIR),C(ICIR),VC(1CIR),DVC(ICIR),CP  B  33 

READ  (1,113)  RD,LD 

READ  (1,140)  STEPFA,OJMAXM,DIELK,DENSPL,DENSMI ,DENS2,REC,EGAP  B  34 

READ  (1,160)  (NXPRNT(I),I=1,6)  B  35 

READ  (1,160)  M, PRTFRQ,N3S,MOOP,MOOFCH, EXTRAP, MS, IPOINT  B  36 

IF  (N3S.GT.0)  GO  TO  20  B  37 

WRITE(6,400)  B  38A 

CHECK  1=1  B  38B 

N3S=1  B  38C 

20  READ  (1,170)  (TPRNTI ( I ) ,TINC( I ) ,NINC( I ) , l=1,N3S)  B  39 

MP1=M+1  B  40 

READ  (1,180)  (NP(I),I=1,MP1)  B  41 

READ  (1,180)  (NM( l),l=1,MP1)  B  42 

GO  TO  27  B  43 

25  IF(ENDING.EQ,1.0R,ENDING,EQ.3)  READ{1,130)USTAT,DU,V0,0MEGA,T1,V  B  43A 

IF(EN0ING.EQ.2.0R.EN0ING.EQ.3)  READ(1,130)  (R( ICIR) ,L( ICIR) ,C( ICIR  B  43C 

*),VC( ICIR),OVC( ICIR),CP,ICIR=2,3)  B  43D 

27  WRITE(6,210)  B  43E 

WRITE  (6,230)  M0DA,A1 ,A2,A3,A4,A5  B  44 

WRITE  (6,220)  B  45 

WRITE  (6,230)  M00B,B1 ,B2,B3,B4,B5  B  46 

WRITE  (6,360)  B  47 

WRITE  (6,230)  M0DD,D1 ,02,D3 ,D4,05  B  48 

WRITE  (6,225)  B  48A 

WRITE  (6,230)  MOOT, CAT, CBT,PWRM,PWRP, TEMP  B  48B 

WRITE  (6,235)  B  49A 

WRITE  (6,190)  DNSTY,SPCHT,THCND  B  49B 

WRITE  (6,370)  B  49 

WRITE  (6,230)  M0DVP,MUP,C1 ,C2,C3,C4  B  50 

WRITE  (6,375)  B  50A 

WRITE  (6,230)  MODVNP,MUPO,AVP,NOP,VPiM  B  50B 

WRITE  (6,380)  B  51 

WRITE  (6,230)  M0DVM,MUM,C5,C6,C7,C8  B  52 

WRITE  (6,385)  B  52A 

WRITE  (6,230)  MODVNM,MUMO,AVM,NOM,VM1M  B  52B 

WRITE  (6,390)  B  53 

WRITE  (6,190)  C9,C10,C15,C20  B  54 

WRITE  (6,240)  B  55 

WRITE  (6,200)  GAMSEC, TAUSEC, GAMMA, GAMI ,JMO,JPO  B  56 

WRITE  (6.250)  B  57 

WRITE  (6,200)  USTAT,DU,V0,0MEGA,T1,V  B  58 

WRITE  (6,520)  B  58A 

WRITE  (6,510)  V1,0MEGA1,PHI 1,V2,0MEGA2,PHI2, lEFCY  B  58B 

WRITE  (6,260)  B  59 

WRITE  (6,200)  P,D,S,T,DVDT,JMAX  B  60 

WRITE  (6,270)  CP  B  61 

DO  30  ICIR=1,3  B  62 
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30  WRITE  (6,280)  ICIR,R( ICIR) ,L( ICIR) ,C( ICIR) ,VC( ICIR) ,DVC( ICIR)  B  63 

WRITE  (6,213)  RD,LD 

WRITE  (6,300)  B  64 

WRITE  (6,290)  STEPFA, DJMAXM, DIELK, DENSPL, DENSMI ,DENS2,REC,EGAP  B  65 

WRITE  (6,310) 

WRITE  (6,320)  (NXPRNT( I ) , 1=1 ,6) 

WRITE  (6,330)  B  66 

WRITE  (6,340)  M,PRTFRQ,N3S,MODP,MOOFCH,EXTRAP,MS, IPOINT  B  67 

IF(ENDING.LE.3)  GO  TO  60  B  69A 

IF  (M0DFGH.LE.1)  GO  TO  40  B  70 

READ  (1,180)  (DN( l),l=1,MP1)  B  71 

GO  TO  52  B  72 

40  DO  50  1=1, MP1  B  73 

50  0N(l)=0.0  B  74 

52  !F  (M0DT.LE.1)  GO  TO  55  B  74A 

READ  (1,180)  (TEM(I),I=1,MP1)  B  74B 

GO  TO  60  B  74G 

55  MTEM=TEMP  B  74D 

DO  57  1=1, MP1  B  74E 

57  TEM(I)=TEMP  B  74F 

G  GALGULATE  AND  PRINT  TPRINT  LIST  B  75 

60  IF  (M0DP.EQ.1)  GO  TO  70  B  76 

IF  (M0DP.EQ.2)  GO  TO  90  B  77 

WRITE(6,410)  B  78A 

GHEGK1=1  B  78B 

RETURN  B  780 

G  TPRINT  LIST  FOR  PRINTING  AGGORDING  TO  TIME  8  79 

70  IJ=0  B  80 

DO  80  J=1,N3S  B  81 

NINK=NING(J)  B  82 

DO  80  l=1,NINK  B  83 

IJ=!J+1  B  84 

Xl  =  l-1  B  85 

80  TPRINT( IJ)=TPRNTI (J)+XI*TING(J)  B  86 

GO  TO  110  B  87 

G  TPRINT  LIST  FOR  PRINTING  AGGORDING  TO  GURRENT  B  88 

90  TPRINT(1)=0.0  B  89 

IJ=1  B  90 

DO  100  J=1,N3S  B  91 

NINK=NING(J)  B  92 

DO  100  l=1,NINK  B  93 

IJ=IJ+1  B  94 

100  TPRINT( IJ)=TPRNTI (J)*TING(J)*»( 1-1)  B  95 

110  WRITE  (6,350)  IJ, (TPR INT( I ) , 1=1 , I J)  B  96 

RETURN  B  97 

G  B  98 

113  FORMAT  (2E13.7) 

120  FORMAT  (72H1 IDENTIF ICATION  CARD  HEADING  EACH  RUN  B  99 

1  )  B  100 

130  FORMAT  (6E13.7)  B  101 

140  FORMAT  (8E10.3)  B  102 

150  FORMAT  ( 1 1 ,E13.7,4E14.7)  B  103 

160  FORMAT  (2013)  B  104 
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170  FORMAT  (2(2E14.8, 1 3 ,5X) )  B  105 

180  FORMAT  (5E14.8)  B  106 

190  FORMAT  (20X,1P,5E20.7)  B  107 

200  FORMAT  (1P,6E20.6)  B  108 

210  FORMAT  (/16X,4HM0DA, 18X,2HA1 , 18X,2HA2, 18X,2HA3 , 18X,2HA4, 18X,2HA5)  B  109 

213  FORMAT  (11X,9HDI0DE  LEG, IP ,2E20.6) 

220  FORMAT  ( 16X,4HM0DB, 18X,2HB1 , 18X,2HB2, 18X,2HB3, 18X,2HB4, 18X,2HB5)  B  110 
225  FORMAT  ( 16X,4HM0DT, 17X,3HCAT, 17X,3HCBT, 16X,4HPWRM, 16X,4HPWRP , 14X,6  B  110A 
1HTEMP  K)  B  110B 

230  FORMAT  ( I20,1P,5E20.7)  B  111 

235  FORMAT  (33X,7HDENSITY,13X,7HSPEC  HT, 12X,8HTHRM  CND)  B  111A 

240  FORMAT  ( 14X,6HGAMSEC, 14X,6HTAUSEC, 15X,5HGAMMA, 16X,4HGAMI , 17X,3HJM0  B  112 
1,17X,3HJP0)  B  112A 

250  FORMAT  ( 15X,5HUSTAT, 18X,2HDU, 18X,2HV0, 15X,5H0MEGA, 18X,2HT1 , 19X  B  113 

I. IHV) 

260  FORMAT  ( 19X, 1HP , 19X, 1HD, 19X, 1HS, 19X, 1HT, 16X,4HDVDT, 16X,4HJMAX)  B  114 

270  FORMAT  ( 15X,27HCIRCUIT  PARAMETERS  WITH  CP=, IP ,E20.6/16X,4H ICIR, 19  B  115 

1X,1HR,19X,1HL,19X,1HC,18X,2HVC,17X,3HDVC)  B  116 

280  FORMAT  ( I20,1P,5E20.6)  B  117 

290  FORMAT  (1P,8E15.4)  B  118 

300  FORMAT  (9X,6HSTEPFA,9X,6HDJMAXM,10X,5H0IELK,9X,6HDENSPL,9X,6HDENSM  B  119 

I I , 10X,5HDENS2,12X,3HREC,11X,4HEGAP)  B  120 

310  FORMAT  (4X,6HNXPRNT)  B  121 

320  FORMAT  (2013)  B  122 

330  FORMAT  ( 11X,1HM,6X,6HPRTFRq,9X,3HN3S,8X,4HM0DP,6X,6HM0DFCH,6X,  B  123 

16H.EXTRAP,10X,2HMS,6X,6HIPOINT)  B  124 

340  FORMAT  (9112)  B  125 

350  FORMAT  (12H  TPRINT  LIST, 1 10, 9H  VALUES/( 1P,8E15.5) )  B  126 

360  FORMAT  ( 16X,4HM00D, 18X,2HD1 , 18X,2HD2, 18X,2HD3, 18X,2HD4, 18X,2HD5)  B  127 

370  FORMAT  ( 15X,5HMODVP , 17X,3HMUP , 18X,2HC1 , 18X,2HC2, 18X,2HC3 , 18X,  B  128 

12HC4) 

375  FORMAT  ( 14X,6HM0DVNP, 16X,4HMUP0,17X,3HAVP, 17X,3HN0P , 16X,4HVP1M)  B  128A 
380  FORMAT  ( 15X,5HMODVM, 17X,3HMUM, 18X,2HC5, 18X,2HC6, 18X,2HC7, 18X,  B  129 

12HC8) 

385  FORMAT  ( 14X,6HM0DVNM, 16X,4HMUM0, 17X,3HAVM, 17X,3HN0M, 16X,4HVM1M)  B  129A 
390  FORMAT  (38X,2HC9,17X,3HC1O,17X,3HC15,17X,3HC20)  B  130 

400  FORMAT  (75H  N3S  (NO.  OF  TPRINT  TRIPLES)  IS  NOT  POSITIVE  AS  REQUIR  B  131 

1ED.  END  OF  PROBLEM.)  B  132 

410  FORMAT  (57H  MODP  IS  NOT  EQUAL  TO  1  OR  2  AS  REQUIRED.  END  OF  PROBL  B  133 

1EM)  B  134 

500  F0RMAT(6E13.7, 12)  B  135 

510  FORMAT(1P,6E20.6,I10)  B  136 

520  F0RMAT(18X,2HV1,14X,6H0MEGA1,16X,4HPHI1,18X,2HV2,14X,6H0MEGA2,16X,  B  137 

*4HPHI2,5X,5HIF0UR)  B  138 

END  B  139- 

SUBROUTINE  INCON(CIRCT, ENDING)  C  1 

REAL  MUM,MUP,NP,NM,JP,JC,JPO,JMO,JM  , I INT, UNTO, IEXT,JINT,JMAX, 123  C  2 

1, ICUR, ICURO,L,JPDIF,JMDIF,JDISPL,MJP,MJM,MJPDF,MJMDF,MJDSP,MTEM,  C  3 

2MUM0,MUP0,N0M,N0P  G  4 

INTEGER  CIRCT(10),ENDING  C  5 

DIMENSION  0NP(101),ONM(101)  C  6 

COMMON  /ABC/  M0DA,M0DB,M0DD,M0DVP,M0DVM,M0DVNP,M0DVNM,M0DT,A1 ,A2,  C  7 


1A3,A4,A5,A6,B1,B2,B3,B4,B5,B6,D1,D2,D3,D4,D5,D6,C1,C2,C3,C4,C5,C6,  C  8 
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5 


6 


C 


C 


C 

10 

20 

C 


C 

30 

C 

40 

50 


2C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,C18,C19,C20,MUP,MUM,MUP0, 
3MUM0 , AVP , AVM , NOP , NOM , VP  1 M , VM 1 M , CAT , CBT , PWRP , PWRM , TEMP 
COMMON  /C I RT/  UST AT , U , V , DU , VO , V 1 , V2 , OMEGA , OMEG A 1 , 0MEGA2 , P H 11 , P H 1 2 , 
1RD,LD,R(3),L(3),C(3),CP,T1, lEFCY.VDD 

COMMON  /DSTR/  NP(201 ) ,NM(201 ) ,DN(201 ) ,TEM(201 ) ,E{201 ) ,DE|201 ) ,VSUM 
1(201),VP(201),VM(201).JP(201),JM(201),JPDIF(201),JMDIF(201),JDISPL 
2(201) ,RCMBR(201) ,ALPHA(201) ,BETA(201),TPRINT(201) 

COMMON  /FUND/  M,MP1 ,N3S,PRTFRQ,NXPRNT(6) ,M0DP,M0DFCH,EXTRAP,MS, I J, 
1 1 PO I  NT , FM , STEPF A , DX , DT , T ,D , S , P , EL , EPSO ,D I  ELK . DNSTY .SPCHT , THCND ,DEN 
2SPI.,DENSMI ,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI , GAMMA, TAUSEC 
COMMON  /M JV/  JPO , JMO , JMAX ,0 JMAXM , JC , J I  NT , I  I  NT , I  I NTO , I  EXT , 1 23 , D I DT , 
1DVDT,DVDTO,VC(3),VCO(3),DVC(3),DVCO(3),ICUR(3),  ICUR0(3) ,MJP,MJM, 
2MJPDF,MJMDF,MJDSP,MTEM,VPMAX,VMMAX,CST(10,3),DMTEM,DDIDT,DIDT0 
IF(ENDING.GT.3)  GO  TO  6 
DO  5  1=1, MP1 
NP{ I )=ONP( I ) 

NM( 1 )=ONM( I ) 

T=OT 

DVOT=ODVDT 

VC(1)=0VC1 

DVC(1)=0DVC1 

GO  TO  55 

EL=1.6022E-19 

EPS0=8.854E-14 

FM=M 

0X=0/FM 

A6=ABD ( MOD A , A 1 , A2 , A3 , A5 ) 

B6=ABD(M0DB,B1,B2,B3,B5) 

D6=ABD(M0DD,D1,D2,D3,D5) 

IF  (MODVP.EQ.O)  C9=SQRT(C3)**3-C3**2/C2*MUP*(1.-C1*C3) 

IF  (MODVM.EQ.O)  C10=SQRT(C7)**3-C7**2/C6*MUM*( 1 ,-C5*C7) 

WRITE  (6,270) 

WRITE  (6,280)  A6,B6,D6,C9,C10 

DETERMINE  TYPE  OF  SERIES  CIRCUIT 

DO  50  ICIR=1,3 

IF  (L( ICIR).Eq.O.)  GO  TO  20 

IF  (C( ICIR) .GE. 1 .E20)  GO  TO  10 

LC  OR  LCR  CIRCUIT 

CIRCT( ICIR)=3 

GO  TO  50 

LR  OR  L  CIRCUIT 

CIRCT( ICIR)=4 

GO  TO  50 

IF  (C( ICIR) .LT.1.E20)  00  TO  30 

IF  (R( ICIR) .NE.O.)  GO  TO  40 

OPEN  CIRCUIT 

CIRCT( ICIR)=1 

GO  TO  50 

RC  CIRCUIT 

CIRCT( ICIR)=2 

GO  TO  50 

R  CIRCUIT 

CIRCT( ICIR)=5 

CONTINUE 


C  9 
C  10 
C  11 
C  12 
C  13 
C  14 
C  15 
C  16 
C  17 
C  18 
C  19 
C  20 
C  21 
C  23 
C  23A 
C  23B 
C  23C 
C  230 
C  23P 
C  23Q 
C  23R 
C  23U 
C  23V 
C  24 
C  25 
C  26 
C  27 
C  28 
C  29 
C  30 
C  31 
C  32 
C  33 
C  34 
C  35 
C  36 
C  37 
C  38 
C  39 
C  40 
C  41 
C  42 
C  43 
C  44 
C  45 
C  46 
C  47 
C  48 
C  49 
C  50 
C  51 
C  52 
C  53 
C  54 
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C  roMPUTE  CIRCUIT  CONSTANTS  C  55 

55  DO  80  ICIR=1,3  C  56 

IF  (R(  ICIR)<tC(ICIR).Eq.O.)  GO  TO  60  C  57 

CST(1,ICIR)=1./(R(ICIR)*C(ICIR))  C  58 

60  IF  (L(ICIR)»C(ICIR).Eq.O.)  GO  10  70  C  59 

CST(2,ICIR)=1./(L(ICIR)<‘C(ICIR))  C  60 

70  IF  (L(ICIR).Eq.O.)  GO  TO  80  C  61 

CST(3, ICIR)=R( ICIR)/L( ICIR)  C  62 

80  CONTINUE  C  63 

CST(4,1)=S»EPS0»DIELK/D  C  64 

CST(5,1)=1./(CP+CST(4,1))  C  65 

CST(6,1)=DX/(2.»EPS0*DIELK)  C  66 

CST(7,1)=EL»CST(6,1)  C  67 

IF(ENDING.LE.3)  GO  TO  90  C  67A 

IF  (MODVP.Eq.1)  CALL  TESTSq  (C1,C2,C3,C9,C15,MUP,C11,C12,C13,C14)  C  68 

IF  (MOOVM.Eq.1)  CALL  TESTSq  {C5,C6,C7,C10,C20,MUM,C16,C17,C18,C19)  C  69 

WRITE  (6,290)  C11 ,C12,C13,C14,C16,C17,C18,C19  C  70 

C  PRINT  INITIAL  DENSITIES  AND  TEMPERATURE  C  71 

90  write  (6,250)  T,V  C  83 

DO  110  1=1, MP1  C  84 

J=l-1  C  85 

110  write  (6,260)  J,NP(I),NM(I),DN(I),TEM(I)  C  86 

C  COMPUTE  U  C  105 

U=USTAT  C  106 

IF  (T.LT.TI)  GO  TO  140  C  107 

IF  (lEFCY.Eq.O)  GO  TO  130  C  108 

U=U+DU+V0*00S(0MEGA*T)+V1*C0S(0MEGA1*T+PHI  1)+V2<»C0S(0MEGA2n+PHI2)  C  109 

GO  TO  140  C  110 

130  U=U+OU+VO*SIN(OMEGA*T)  C  111 

140  DO  150  1=1, MP1  C  112 

150  J0ISPL(  l)=DIELK<»EPSO*DVDT/D  C113 

C  STORE  INITIAL  CONDITIONS  C  115 

DO  220  ICIR=1,3  C  117 

0VC0( ICIR)=DVC( ICIR)  C118 

VCO( ICIR)=VC( ICIR)  C  119 

VS=V  C  120 

IF  ( ICIR.EQ.1)  VS=U-V  C  121 

ICIRT=CIRCT( ICIR)  C  122 

GO  TO  (190,200,200,190,210),  ICIRT  C  123 

190  ICUR( ICIR)=VC( ICIR)  C  124 

GO  TO  220  C  125 

200  ICUR( ICIR)=C( ICIR)*0VC( ICIR)  C  126 

GO  TO  220  C  127 

210  ICUR( ICIR)=VS/R( ICIR)  C128 

220  ICURO( ICIR)=ICUR( ICIR)  C129 

IEXT=ICUR(1)  C  131 

DIDT=0  C  133 

DDIDT=0 

I23=ICUR(2)-+ICUR(3)  C  134 

IF(ENDING.LE,3)  RETURN  C  134A 

DO  230  1=1, MP1  C  142 

ONP( I )=NP( : )  C  143 

230  ONM(l)=NM(l)  C  144 

0T=T  C  156 
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ODVDT=DVDT  C  157 

0VC1=VC(1)  C  158 

0DVC1=DVC(1)  C  159 

RETURN  C  162 

250  FORMAT  (3H1T=,1PE16.6,5H  V=,E16 .6//4X, 1HM, 18X,2HN+, 18X,2HN-, 18X,  C  163 

12HDN,14X,6HTEMP  K)  C  164 

260  FORMAT  ( I5,1P,4E20.6)  C  165 

270  FORMAT  (/38X,2HA6,18X,2HB6,18X,2HD6,18X,2HC9,17X,3HC10)  C  166 

280  FORMAT  (20X,1P,5E20.6)  C  167 

290  FORMAT  ( 37X,3HC1 1 , 17X,3HC12, 17X,3HC13 , 17X,3HC14/20X, IP ,4E20.7/37X,  C  168 

13HC16,17X,3HC17,17X,3HC18,17X,3HC19/20X,4E20.7)  C  169 

END  C  170- 

SUBROUTINE  EXTCIR  (ISET,CIRCT)  D  1 

REAL  JC,I INT,JINT,I INTO,  MJDSP, IEXT,L, ICUR, ICURO, 123, UMAX, JMO,JPO  D  2 

COMMON  /GIRT/  USTAT,U,V,DU,V0,V1 ,V2,0MEGA,0MEGA1 ,0MEGA2,PH 1 1 ,PHI2,  D  3 

1RD,LD,R(3),L(3),C(3),CP,T1,IEFCY,VDD  D  4 

COMMON  /FUND/  M,MP1 ,N3S,PRTFRQ,NXPRNT(6) ,MODP,MODFCH, EXTRAP, MS, IJ,  D  5 

1IP0INT,FM,STEPFA,DX,DT,T,D,S,P,EL,EPS0,DIELK,DNSTY,SPCHT,THCND,DEN  D  6 

2SPL,DENSMI,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI, GAMMA, TAUSEC  D  7 

COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT,l INT,I INT0,IEXT,I23,DIDT,  D  8 

1DVDT,DVDTO,VC(3),VCO(3),OVC(3),DVCO(3),ICUR(3),  ICUR0(3) ,MJP,MJM,  D  9 

2MJPDF ,MJMDF ,MJDSP ,MTEM, VPMAX, VMMAX ,CST( 10,3) ,DMTEM,DD I DT,D I DTO  D  10 

REAL  MJM,MJP,MJMDF,MJPDF,MTEM  D  11 

INTEGER  CIRCT(IO)  D  12 

IF  ( ISET.EQ.2)  GO  TO  20  D  13 

C  UPDATE  QUANTITIES  0  14 

IINTO=IINT  D  15 

DIDTO=DIDT 

DO  10  ICIR=1,3  D  16 

DVC0( ICIR)=DVC( ICIR)  D  17 

VCO( ICIR)=VC( ICIR)  0  18 

10  ICURO( ICIR)=ICUR( ICIR)  D  19 

C  COMPUTE  CURRENT  IN  SERIES  CIRCUITS  D  25 

20  DO  90  ICIR  =1,3  0  26 

VS=V  D  27 

IF  (ICIR.EQ.I)  VS=U-V  D  28 

ICIRT=CIRCT( ICIR)  D  29 

GO  TO  (40,50,60,70,80),  ICIRT  D  30 

40  ICUR{ICIR)=0  D  31 

GO  TO  90  D  32 

50  DVC{ ICIR)=CST(1,ICIR)*(VS-VC( ICIR))  D  33 

ICUR( ICIR)=C( ICIR)*DVC( ICIR)  D  34 

GO  TO  90  D  35 

60  DVC2=CST(2, ICIR)*(VS-VC( ICIR) )-CST(3, ICIR)*DVC( ICIR)  D  36 

0VC( ICIR)=0VC0{ ICIR)+DT*DVC2  D  37 

ICUR( ICIR)=C( ICIR)*DVC(  ICIR)  D  38 

GO  TO  90  D  39 

70  DI=(1./L( ICIR))*VS-CST(3,ICIR)*ICUR0( ICIR)  D  40 

ICUR( ICIR)=ICUR0( ICIR)+DT*DI  D  41 

GO  TO  90  0  42 

80  ICUR( ICIR)=VS/R( ICIR)  D  43 

90  VC(  ICIR)=VC0(  ICIR)+DT<»DVC(  ICIR)  D  44 

C  COMPUTE  VOLTAGE  AND  CURRENT  OF  DEVICE  D  45 
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c 

20 


50 


C 

60 

70 

80 

C 

90 

C 


DVDT=( ICUR(1)-ICUR(2)-ICUR(3)-I INT)/CP 

V=V+OT«DVDT 

IEXT=ICUR(1) 

RETURN 

END 

SUBROUT  I NE  TEST( M I NUS , I NSTEP .CHECK ,CHECK2 , 1 SET ,MSO) 

REAL  NMSAV.NPSAV.JC.MJDSP,  UMAX,  I INT.JINT, UNTO, 123, ICUR, ICURO, 
<*IEXT,NP,NM.JP,JM,JCSAV,JDISPL,JMDIF,JPDIF,MJP,MJM,MJPDF,MJMDF 
INTECER  CHECK, CHECK2 

COMMON  /DSTR/  NP(201 ) ,NM(201 ) ,DN{201 ) ,TEM(201 ) ,E(201 ) ,DE{201 ) ,VSUM 
1(201),VP(201),VM(201),JP(201),JM(201),JPDIF(201),JMDIF(201),JDISPL 
2(201) ,RCMBR(201) ,ALPHA(201) ,BETA(201),TPRINT(201) 

COMMON  /FUN0/M,MP1,N3S,PRTFRQ,NXPRNT{6),M0DP,M0DFCH, EXTRAP , MS , IJ , 
1IPOINT,FM,STEPFA,DX,OT,T,D,S.P,EL.EPSO,DIELK,DNSTY,SPCHT,THCND,DEN 
2SPL,DENSMI ,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI .GAMMA, TAUSEC 
COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT, I INT, I  INTO, lEXT, I23,DIDT, 
1DVDT,DVDTO,VC(3),VCO(3),OVC(3),OVCO(3),ICUR{3)  ,  ICURO( 3 ) ,MJP ,MJM, 
2MJPDF,MJMDF,MJDSP,MTEM,VPMAX,VMMAX,CST(10,3),DMTEM,DDIDT,DIDTO 
COMMON  /S AV/  NPS AV ( 20 1 ) , NMSAV (201), ESAV (201), VXSAV ( 20 1 ) , VPSAV ( 20 1 ) 
1 , VMSAV ( 20 1 ) , JPSAV ( 20 1 ) , JMSAV (201), VSAV, TSAV, JCSAV ,GAMSAV , TMPSAV 
REAL  JMSAV, JPSAV, JMO,JPO,MTEM 
IF  ( INSTEP. NE.O)  GO  TO  60 
STORE  PRESENT  VALUE  OF  VARIABLES 
DO  50  1=1, MP1 
NPSAV( l)=NP( I) 

NMSAV( I )=NM(  1 ) 

VPSAV(I)=VP(I) 

VMSAV(I)=VM(I) 

JPSAV( I )=JP( I ) 

JMSAV( I )=JM( I ) 

VXSAV ( I )=VSUM( I) 

ESAV( I )=E(  I ) 

TSAV=T 

VSAV=V 

GAMSAV=GAMMA 

JCSAV=JC 

TMPSAV=MTEM 

ISET=1 

IF  ( CHECK. EQ.O)  MINUS=MSO 
GO  TO  150 

CHECK  FOR  NEGATIVE  VELOCITIES 
DO  70  1=1, MP1 

IF  (VM(I).LT.O..OR.VP(I).LT.O.)  GO  TO  80 

CONTINUE 

GO  TO  90 

CONTINUE 

CHECK=1 

CHECK  FOR  LARGE  CURRENT  OR  TIME 
IF  (M0DP.EQ.2)  GO  TO  100 
IF  (ABS(JC) .LE.JMAX)  GO  TO  110 
CURRENT  DENSITY  TOO  LARGE 
WRITE  (6,160) 

GO  TO  140 


D  46 
D  47 
D  48 
D  53 
D  54- 
E  1 
E  2 
E  3 
E  4 
E  5 
E  6 
E  7 
E  8 
E  .9 
E  10 
E  11 
E  12 
E  13 
E  14 
E  15 
E  17 
E  22 
E  37 
E  38 
E  39 
E  40 
E  41 
E  42 
E  43 
E  44 
E  47 
E  48 
E  49 
E  51 
E  52 
E  52A 
E  52B 
E  53 
E  54 
E  55 
E  56 
E  57 
E  58 
E  59 
E  60 
E  61 
E  62 
E  63 
E  64 
E  65 
E  67 
E  67A 
E  68 
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100  IF  (T.LE.JMAX)  GO  TO  110  E  69 

WRITE  (6,170)  E  70 

GO  TO  140  E  71 

C  CHECK  FOR  TOO  RAPID  CURRENT  CHANGE  E  72 

110  DJC=ABS((JC-JCSAV)/JC)  E  73 

IF  (DJC.LE.DJMAXM)  GO  TO  20  E  73A 

WRITE  (6,180)  E  738 

WRITE  (6,190)  T,JC.JCSAV,DJC  E  73C 

MINUS=MSO  E  74 

CHECK=0  E  75 

CHECK2=CHECK2+1  E  77 

IF  (CHECK2.LE.2)  GO  TO  120  E  78 

GO  TO  140  E  80 

C  RESET  ARRAYS  TO  SAVE  VALUES  AND  REPEAT  INSTEP  E  81 

120  DO  130  1=1, MP1  E  82 

NP( I )=NPSAV( I )  E  83 

NM{ I )=NMSAV( I )  E  84 

VP( I )=VPSAV( I )  E  85 

VM( I )=VMSAV( I )  E  86 

JP( I )=JPSAV( I )  E  87 

JM( I )=JMSAV( I )  E  88 

VSUM( I )=VXSAV( I )  E  91 

130  E(I)=ESAV(I)  E  92 

T=TSAV  E  93 

V=VSAV  E  95 

GAMMA=GAMSAV  E  96 

JC=JCSAV  E  96A 

MTEM=TMPSAV  E  968 

ISET=2  E  97 

GO  TO  150  E  98 

140  ISET=3  E  99 

150  RETURN  E  100 

160  FORMAT  (/29H  CURRENT  DENSITY  EXCEEDS  UMAX)  E  101 

170  FORMAT  (/18H  TIME  EXCEEDS  UMAX)  E  102 

180  FORMAT  (/33H  TIME  STEP  REDUCED, T,JC,JCSAV,DJC)  E  103 

190  FORMAT  (1P,4E20.7)  E  103A 

END  E  104- 

SUBROUTINE  OUTPUT(DVXDT)  F  1 

REAL  MJP,MJM,JDIF,MJDSP,IEXT,I INT,JINT,I23,L,I INTO,JMO,JPO,JC,JP,  F  2 
1JM,JMAX, ICUR, ICURO,MJMDF,MJPDF,JPDIF,JMOIF,JDISPL,NP,NM,MTEM  F  3 

INTEGER  PRTFRQ  F  4 

DOUBLE  PRECISION  TITLE(40) 

DIMENSION  XPR(6),DVXDT(101) 

COMMON  /DSTR/  NP (201 ) ,NM(201 ) ,DN( 201 ) ,TEM(201 ) ,E(201 ) ,DE(201 ) ,VSUM  F  6 
1(201) ,VP(201) ,VM(201) ,JP(201) ,JM(201),JPDIF(201) ,JMDIF(201) ,JDISPL  F  7 
2(201) ,RCI1BR(201),ALPHA(201),BETA(201),TPRINT(201)  F  8 

COMMON  /FUND/  M,MP1 ,N3S, PRTFRQ, NXPRNT(6) ,MOOP,MODFCH, EXTRAP, MS, IJ,  F  9 
1 IPOr;T,FM,STEPFA,DX,DT,T,D,S,P,EL,EPSO,DIELK,DNSTY,SPCHT,THCND,DEN  F  10 
2SPL,DENSMI ,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI , GAMMA, TAUSEC  F  11 

COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT,l INT, I  INTO, lEXT, I23,DIDT,  F  12 

1DV01,DVDTO,VC(3),VCO(3),DVC(3),DVCO(3),ICUR(3),  ICUR0(3) ,MJP,MJM,  F  13 
2MJPDF,MJMDF,MJDSP,MTEM,VPMAX,VMMAX,CST(10,3) ,DMTEM,DDIDT,DIDTO  F  14 

COMMON  /CIRT/  USTAT,U,V,DU,V0,V1 ,V2,0MEGA,0MEGA1 ,0MEGA2,PH 1 1 ,PHI2,  F  15 
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1RD,LD,R(3),L(3),C(3),CP,T1,IEFCY,VNEW  F  16 

DATA  TITLE  /6H  A,4HLPHA,6H  ,4HBETA,6H  D,4H  F  20 

1ELTA,6H  V,4H-(X),6H  V,4H+(X),6H  ALPH,4HA*J-,6H  BET,4HA*  F  21 

2J+,6H  N,4HET  Q,6HREC0M  ,4HRATE,6H  J+(D,4HIFF) ,6H  J-(D,4HIFF)  F  22 

3,6H  J  0,4HISPL,6H  J  T,4HOTAL,6H  ,4HV(X),6H  DV{X,4H)/DT,6  F  23 

4H  N,4H+(X),6H  N,4H-(X),6K  J,4H+(X),6H  J,4H-(X),6H  F  24 

5  TE,4HMP  K/  F  24A 

UNIVD(EA)=UNIV(P.M0DD,EA,D1,D2,D3,D4,D5,D6)  F  25 

C  COMPUTE  MEAN  VALUES  F  26 

BJPDX=(BETA(1)*JP(l)+BETA(MP1)*JP{MP1))/2.  F  29 

AJMDX=(ALPHA{1)*JM(1)+ALPHA(MP1)*JM(MP1))/2.  F  30 

BETDX=(BETA(1)+BETA(MP1))/2.  F  31 

ALPDX=(ALPHA(1)+ALPHA(MP1))/2,  F  32 

RCMBAV=(RCMBR(1)+RCMBR(MP1))/2.  F  33 

DO  10  1=2, M  F  34 

RCMBAV=RCMBAV+RCMBR( I)  F  35 

ALPOX=ALPDX+ALPHA( I )  F  36 

BETDX=BETDX+B£TA( I)  F  37 

AJM0X=AJM0X+ALPHA( I ) « JM( I )  F  38 

10  BJPOX=BJPOX+BETA( l)*JP( I)  F  39 

ALPDX=ALPDX*DX  F  42 

BETOX=BETOX»DX  F  43 

AJMDX=AJMDX/FM  F  44 

RCMBAV=RCMBAV/FM  F  45 

BJPDX=BJPDX/FM  F  48 

JD I F=MJPOF+MJMOF  F  49 

I1=2**NXPRNT(1)-1  F  53 

I2=2<*NXPRNT(2)-1  F  54 

I3=2*NXPRNT(3)-1  F  55 

I4=2*NXPRNT(4)-1  F  56 

(5=2«NXPRNT(5)-1  F  57 

I6=2»NXPRNT(6)-1  F  57A 

WRITE  (6,240)  T,VNEW,U,MJP,MJM,  JC  ,V,OT,VPMAX,VMMAX, lEXT, I INT, J  F  58 

1 1  NT , DVDT , ALPDX , A JMOX , BETDX , B JPDX ,RCMB AV , 1 23 , D I DT , H JPDF , M JMDF , JD I F ,  F  59 

2MJDSP,MTEM,DMTEM,DENS2T,TITLE( II  F  60 

3) ,TITLE(I1+1),TITLE( I2),TITLE( 12+1) ,TITLE( 1 3) ,TITLE( 1 3+1 ) ,TITLE( 14  F  61 

4) ,TITLE{I4+1),TITLC; I5),TITLE( 15+1) ,TITLE( 16) ,TITLE( 16+1)  F  61A 

VSUM(1)=0.  F  62 

XPRT=0.  F  63 

XMPLY=PRTFRQ  F  64 

DO  230  l=1,MP1,PRTFRQ  F  65 

DO  220  J=1,6  F  66 

NJ=NXPRNT(J)  F  67 

GO  TO  (30,40,50,60,70,80,90,100,110,120,130,140,150,160,170,180,19  F  68 

10,200,210,215),  NJ  F  69 

30  XPRINT=ALPHA( I )  F  70 

GO  TO  220  F  71 

40  XPRINT=BETA(I)  F  72 

GO  TO  220  F  73 

50  XPRINT=UNIVD(E(I))  F  74 

GO  TO  220  F  75 

60  XPRINT=VM(I)  F  76 

GO  TO  220  F  77 
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70 

XPRINT=VP( 1 ) 

F 

78 

GO  TO  220 

F 

79 

80 

XPRINT=ALPHA(I)*JM(I) 

F 

80 

GO  TO  220 

F 

81 

90 

XPRINT=BETA(I)*JP(I) 

F 

82 

GO  TO  220 

F 

83 

100 

XPRINT=NP( 1 )-NM( 1 )+DN(  1 ) 

F 

84 

GO  TO  220 

F 

85 

110 

XPRINT=RCMBR( 1) 

F 

88 

GO  TO  220 

F 

89 

120 

XPRINT=JPDIF( 1 ) 

F 

90 

GO  TO  220 

F 

91 

130 

XPRINT=JMDIF( 1 ) 

F 

92 

GO  TO  220 

F 

93 

140 

XPRINT=JDISPL( 1 ) 

F 

94 

GO  TO  220 

F 

95 

150 

XPR 1 NT= JP ( 1 ) + JM ( 1 ) + JO  1 SPL ( 1 ) + JPD 1 F { 1 ) +JMD 1 F ( 1 ) 

F 

96 

GO  TO  220 

F 

97 

160 

XPRINT=VSUM( 1 ) 

F 

98 

GO  TO  220 

F 

99 

170 

XPRINT=OVXDT( 1 ) 

F 

100 

GO  TO  220 

F 

101 

180 

XPRINT=NP( 1 ) 

F 

102 

GO  TO  220 

F 

103 

190 

XPRINT=NM( 1 ) 

F 

104 

GO  TO  220 

F 

105 

200 

XPRINT=JP( 1 ) 

F 

.106 

GO  TO  220 

F 

107 

210 

XPRINT=JM( 1 ) 

F 

108 

GO  TO  220 

F 

108A 

215 

XPRINT=TEM(I) 

F 

108B 

220 

XPR( J)=XPRINT 

F 

109 

WRITE  (6,250)  XPRT,E( 1 ) , (XPR( J) , J=1 ,6) 

F 

11.0 

230 

XPRT=XPRT+XMPLY*DX 

F 

111 

RETURN 

F 

112 

C 

F 

113 

240 

FORMAT  ( 1 H 1 , 1 4X , 1 HT , 12X ,4HVNEW , 15X , 1 HU, 1 IX , 5H J+AVG , 1 IX , 5H J- AVG , 1 1 

F 

114 

1X,5HJ  AVG,13X,1HV/1P7El6.6/14X,2HDT,1tX,5HV+MAX,11X,5HV-MAX,11X,5 

F 

115 

2HI  EXT,11X,5HI  INT,11X,5HJ  INT.11X,5HDV/DT/7E16.6/6X,10H0»ALPHAAVG 

F 

116 

3 ,9X,7HA*J-AVG,7X,9HD»BETAAVG, 9X,7HB»J+AVG,5X, 1 1HREC0MRATAVG, 10X, 

F 

117 

46H  I23,11X,5H0I/DT/1P,7E16.6/11X,5HJ+0IF,11X,5HJ-DIF,11X, 

F 

118 

55HJ  DIF,9X,7HJ  DISPL,6X,10HAVE  TEMP  K,8X,8HTEM  RATE,10X, 
66HDENS2T/1P,7E16.6/9X,1HX,13X,4HE(X) 

7,6(7X,A6,A4)) 

F 

119 

250 

FORMAT  (2X,1PE10.3,1P,7E17.5) 

F 

121 

END 

F 

122- 

FUNCTION  ABD  (M0D,H1,H2,H3,H4) 

H 

1 

C 

A6  CALLS  ABD(M0DA,A1,A2,A3,A5) 

H 

2 

C 

B6  CALLS  ABD(M0DB,B1,B2,B3,B5) 

H 

3 

C 

D6  CALLS  ABD(M0DD,D1, 02,03,05) 

H 

4 

ABC=0 . 

H 

5 

IF  (H4.EQ.0. .OR.MOD.EQ.O)  GO  TO  10 

H 

6 

IF  (M0D.EQ.2)  H4=SQRT(H4) 

H 

7 

ABC=-H4* { ALOG ( H 1 /H3 ) -H2/H4 ) 

H 

8 
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10  ABD=ABC  H  9 

RETURN  H  10 

END  H  11 

FUNCTION  UNIV  (XP,MOD,EA,G1 ,G2,G3.G4,G5,G6)  I  1 

C  UNIVA(EA)  =  UNIV(XP,MODA,EA,A1,A2,A3.A4,A5,A6)  I  2 

C  UNIVB(EA)  =  UNIV(XP.MODB.EA,B1,B2,B3,B4,B5,B6)  I  3 

C  UNIVD(EA)  =  UNIV(XP,M0D0,EA,D1,D2,D3,D4,D5,D6)  I  4 

VAR1(F1,F2)=XP*F1*EXP(-F2*(XP/EB)**XP0N)  I  5 

VAR=0.  I  6 

IF  (MOD.EQ.O)  GO  TO  70  17 

EB=ABS(EA)  I  8 

EAP=EB/XP  I  9 

IF  (EAP.LE.G2/100.)  GO  TO  70  I  10 

XP0N=1.  I  11 

GO  TO  (10,40,50,60),  MOD  I  12 

10  IF  (EAP.GT.G5)  GO  TO  30  I  13 

20  VAR=VAR1(G1,G2)  I  14 

GO  TO  70  I  15 

30  VAR=VAR1{G3,G6)  I  16 

GO  TO  70  I  17 

40  XP0N=.5  I  18 

IF  (EAP.GT.G5)  GO  TO  30  I  19 

GO  TO  20  I  20 

50  VAR=VAR1{G1,G2)+VAR1(G3,G4)  I  21 

GO  TO  70  I  22 

60  XPONs.5  I  23 

GO  TO  50  I  24 

70  UNIV=VAR  I  25 

RETURN  I  26 

END  I  27 

REAL  FUNCTION  MUF(Z,P,G1,G2,G3,G4,G5,MU)  J  1 

C  MUPF(Z)  CALLS  MUF(Z,P,C1 ,C2,C3,C4,C9,MUP)  J  2 

C  MUMF(Z)  CALLS  MUF(Z,P,C5,C6,C7,C8,C10,MUM)  J  3 

REAL  MU  J  4 

Y=ABS(Z)  J  5 

IF  (Y,GT.G3*P.ANO.Y.NE.O)  GO  TO  10  J  6 

MUF=MU»{1,-G1<»Y/P)  J  7 

RETURN  J  9 

10  Y=SQRT(1./Y)  J  10 

MUF=G2*Y*(1.-G5*Y**3)  J  11 

20  IF  (MUF.GT.G4*Y»*2)  MUF=G4<*Y**2  J  12 

RETURN  J  13 

END  J  14 

SUBROUTINE  TESTSQ  (01 ,D2,D3 ,D4,D5,D6,E1 ,E2,E3 ,E4)  K  1 

DIMENSION  A(4,4),  B(4)  K  2 

A(1,1)=1.  K  3 

A(1,2)=D3  K  4 

A(1,3)=D3»D3  K  5 

A(1,4)=D3*D3»D3  K  6 

A(2,1)=1.  K  7 

A(2,2)=D5  K  8 

A(2,3)=D5*D5  K  9 

A(2,4)=D5*D5*D5  K  10 
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A(3,1)=0.  K  11 

A(3,2)=1.  K  12 

A(3,3)=2.»D3  K  13 

A(3,4)=3.»D3*D3  K  14 

A(4,1)=0.  K  15 

A(4,2)=1,  K  16 

A(4,3)=2.*D5  K  17 

A(4,4)=3.*D5*D5  K  18 

B(1)=D6*D3<»{1.-D1*D3)  K  19 

B(2)=D2*(SQRT(D5)-D4/D5)  K  20 

B(3)=D6*(1.-2.»D1«D3)  K  21 

B(4)=D2*(1./(2.*SQRT(D5))+04/(D5*D5))  K  22 

CALL  SIMQ  (A,B,4.KS)  K  23 

E1=B{1)  K  24 

E2=B(2)  K  25 

E3=B(3)  K  26 

E4=B(4)  K  27 

RETURN  K  28 

END  K  29- 

SUBROUTINE  SIMQ(A,B.N,KS)  L  48 

DIMENSION  A(1),  8(1)  L  49 

C  L  50 

C  FORWARD  SOLUTION  L  51 

C  L  52 

TOL=0.0  L  53 

KS=0  L  54 

JJ=-N  L  55 

DO  80  J=1,N  L  56 

JY=J+1  L  57 

JJ=JJ+N+1  L  58 

BIGA=0  L  59 

IT=JJ-J  L  60 

DO  20  l=J,N  L  61 

C  L  62 

C  SEARCH  FOR  MAXIMUM  COEFFICIENT  IN  COLUMN  L  63 

C  L  64 

IJ=IT+I  L  65 

IF  (ABS(BIGA)-ABS(A( IJ) ) )  10,20,20  L  66 

10  BIGA=A(IJ)  L  67 

IMAX=I  L  68 

20  CONTINUE  L  69 

C  L  70 

C  TEST  FOR  PIVOT  LESS  THAN  TOLERANCE  (SINGULAR  MATRIX)  L  71 

C  L  72 

IF  (ABS(BIGA)-TOL)  30,30,40  L  73 

30  KS=1  L  74 

RETURN  L  75 

C  L  76 

C  INTERCHANGE  ROWS  IF  NECESSARY  L  77 

C  L  78 

40  l1=J+N*(J-2)  L  79 

IT=IMAX-J  L  80 
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DO  50  K=J,N  L  81 

I1=I1+N  L  82 

I2=I1+IT  L  83 

SAVE=A(I1)  L  84 

A(n)=A(l2)  L  83 

A(t2)=SAVE  L  86 

C  L  87 

C  DIVIDE  EQUATION  BY  LEADING  COEFFICIENT  L  88 

C  L  89 

SAVE=B(IMAX)  L  90 

50  A( I1)=A( I1)/BIGA  L  91 

B(IMAX)=B(J)  L  92 

B(J)=SAVE/BIGA  L  93 

C  L  94 

C  ELIMINATE  NEXT  VARIABLE  L  95 

C  L  96 

IF  (J-N)  60,90,60  L  97 

60  IQS=N*(J-1)  L  98 

DO  80  IX=JY,N  L  99 

IXJ=iqS+IX  L  100 

IT=J-IX  L  101 

DO  70  JX=JY,N  L  102 

IXJX=N*(JX-1)+IX  L  103 

JJX=IXJX+IT  L  104 

70  A( IXJX)*A( IXJX)-(A( IXJ)*A(JJX))  L  105 

80  B( IX)=B( IX)-(B(J)»A( IXJ))  L106 

C  L  107 

C  BACK  SOLUTION  L  108 

C  L  109 

90  NY=N-1  L  110 

IT=N*N  L  111 

DO  100  J=1,NY  L  112 

IA=IT-J  L  113 

IB=N-J  L  114 

IC=:N  L  115 

DO  100  K=1,J  L  116 

B( IB)=B( IB)-A( IA)*B( 1C)  L117 

IA=IA-N  L  118 

100  IC=IC-1  L  119 

RETURN  L  120 

END  L  121- 

SUBROUTINE  EFFCY  M  1 

REAL  L,JAVG,ICUR,ICURO,IEXT,IINT,l INT0,I23,JC,JINT,JMAX,JM0,JP0,  M  2 

1MJM,MJP,MJM0F,MJPDF,MJ0SP,MTEM  M  2A 

COMMON  /CIRT/  USTAT,U,V,DU,V0,V1 ,V2,0MEGA,0MEGA1 ,0MEGA2,PHI 1 ,PHI2,  M  3 

1RD,LD,R(3),L(3),C(3),CP,T1,IEFCY,VDD  M  4 

COMMON  /EFCY/  AMP, AMP2, PHASE, PJ,PV,PHASE2,PJ2,PV2, DTP, GNEG,RFP,EFF  M  5 

1,AMPVDC,AMPJDC,DCPWR,AMPV,AMPV2  M  6 

COMMON  /FUND/  M,MP1 ,N3S,PRTFRQ,NXPRNT(6) ,MODP,MOOFCH, EXTRAP, MS, IJ,  M  7 

1IP0INT,FM,STEPFA,DX,DT,T,D,S,P, EL, EPSO, DIELK, DNSTY,SPCHT,THCND, DEN  M  8 

2SPL,DENSMI,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI, GAMMA, TAUSEC  M  9 

COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT,l INT,I INT0,IEXT,I23,DIDT,  M  10 

1DVDT,DVDTO,VC(3),VCO(3),DVC(3),DVCO(3),ICUR(3),  ICUR0(3) ,MJP,MJM,  M  11 

2M JPDF , M JMDF , M JDSP , MTEM , VPMAX , VMMAX ,CST (10,3), DMTEM ,0D I DT , D I DTO  M  12 
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COMMON  /PLO/  TIME(401),VOLT(401),CURR(401),JAVG(401),CURR1(401),  M  13 

1CURR2(401),CURR3(401),CLTT(401),TAVR(401)  M  14 

DATA  DPR/57.295780/  M  15 

C  CALCULATE  TIME  POINTS  M  16 

TP=2. *3. 14159/OMEGA  M  17 

TK=TIME( IPOINT)-TP  M  18 

DO  10  I1=1,IP0INT  M  19 

I2=IP0INT-I1+1  M  20 

IF(TIME| I2).GT.TK)  GOTO  10  M  21 

GO  TO  20  M  22 

10  CONTINUE  M  23 

20  CJ=JAVG(I2)+(JAVG(I2+1)-JAVG(I2))*(TIME(I2)-TK)/DTP  M  24 

VJ=VOLT( I2)+(V0LT( I2+1)-V0LT( I2))*(TIME{ I2)-TK)/DTP  M  25 

DTP2=TIME( I2+1)-TK  M  26 

CTO=COS(OMEGA*TK)  M  27 

CT1=C0S(0MEGA*TIME( 12+1))  M  28 

CT2=C0S  ( OMEGA^^T  I  ME  (  I PO I  NT ) )  M  29 

CT02=C0S(2.+0MEGA»TK)  M  30 

CT12=C0S(2.*0MEGA+TIME( 12+1))  M  31 

CT22=C0S(2.<»0MEGA*TIME( IPOINT))  M  32 

STO=SIN(OMEGA»TK)  M  33 

ST1=SIN(0MEGA»TIME( 12+1))  M  34 

ST2=SIN(0MEGA<»TIME(  IPOINT))  M  35 

ST02=SIN(2.*OMEGA*TK)  M  36 

ST12=SIN(2.»0MEGAmME(  12+1))  M  37 

ST22®SIN(2.»0ME0A»TIME( IPOINT))  M  38 

SUMJ1=(JAVG( I2+1)»CT1+JAVG( IP0INT)*CT2)/2.  M  39 

SMJ1*{CJ»CT0+JAVG(l2+1)«CT1)/2.  M  40 

SUMJ2=(JAVG( I2+1)*ST1+JAVG( IP0INT)*ST2)/2.  M  41 

SMJ2={CJ+ST0+JAVG( l2+1)»ST1)/2.  M  42 

SUMJ12=(JAVG( I2+1)»CT12+JAVG( IP0INT)»CT22)/2.  M  43 

SM J 1 2= ( C J+CT02+ JAVG ( 1 2+ 1 ) *CT 1 2 ) /2 .  M  44 

SUMJ22=(JAVG( I2+1)*ST12+JAVG( IP0INT)»ST22i/2.  M  45 

SMJ22=(CJ+ST02+JAVG( l2+1)<*ST12)/2.  M  46 

SUMV1=(V0LT( I2+1)<*CT1+V0LT( IP0INT)»CT2)/2.  M  47 

SMV1=(VJ*CTO+VOLT{ l2+1)*CT1)/2.  M  48 

SUMV2= ( VOLT ( 1 2+ 1 ) +ST 1 +VOLT ( I PO I  NT ) *ST2 ) /2 .  M  49 

SMV2=:(VJ*ST0+V0LT{ l2+1)*ST1)/2.  M  50 

SUMV12=(V0LT( I2+1)*CT12+V0LT( IP0INT)+CT22)/2.  M  51 

SMV12=(VJ*CT02+VOLT( l2+1)*CT12)/2.  M  62 

SUMV22=(V0LT( I2+1)*ST12+V0LT( IP0INT)»ST22)/2.  M  53 

SMV22= ( V J*ST02+VOLT ( 1 2+ 1 ) *ST 1 2 ) /2 .  M  54 

SUMJD=(JAVG( I2+1)+JAVG( IP0INT))/2.  M  55 

SMJD=(CJ+JAVG( 12+1) )/2.  M  56 

SUMVD=(VOLT{ I2+1)+V0LT( IP0INT))/2.  M  57 

SMV0=  (VJ+VOLT{ l2+1))/2.  M  58 

122=12+2  M  59 

IP1=IP0INT-1  M  60 

DO  30  J1«I22,IP1  M  61 

CT=C0S(0MEGA*TIME(J1))  M  62 

ST=SIN(0MEGA*TIME(J1))  M  63 

CT2=C0S(2.+0MEGA*TIME(J1))  M  64 

ST2=SIN(2.«0MEGA*TIME(J1))  M  65 

SUMJ1=JAVG( J1)*C1+SUMJ1  M  66 
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SUMJ2=JAVG( J1)*ST+SUMJ2  M  67 

SUMJ12=JAVG( J1)*CT2+SUMJ12  M  68 

SUMJ22=JAVG{J1)*ST2+SUMJ22  M  69 

SUMV1=V0LT{J1)»CT+SUMV1  M  70 

SUMV2=V0LT(J1)*ST+SUMV2  M  71 

SUMV12=V0LT(J1)»CT2+SUMV12  M  72 

SUMV22=V0LT ( J 1 ) *ST2+SUMV22  M  73 

SUMJD=JAVG(J1)+SUMJD  M  74 

30  SUMVD=V0LT(J1)+SUMVD  M  75 

SUMJ1=2,*(SUMJ1*DT?'+SMJ1*DTP2)/TP  M  76 

SUMJ2=2.*(SUMJ2»DTP+SMJ2»DTP2)/TP  M  77 

SUMJ12=2.*(SUMJ12<fDTP+SMJ12*DTP2)/TP  M  78 

SUMJ22=2.*(SUMJ22*DTP+SMJ22*DTP2)/TP  M  79 

SUMV1=2.*(SUMV1*DTP+SMV1*DTP2)/TP  M  80 

SUMV2=2.*(SUMV2*0TP+SMV2*DTP2)/TP  M  81 

SUMV12=2.*(SUMV12»DTP+SMV12*DTP2)/TP  M  82 

SUMV22=2.<*(SUMV22*DTP+SMV22»DTP2)/TP  M  83 

AMPJ0C=(SUMJD*DTP+SMJD»DTP2)/TP  M  84 

AMPVDC=(SUMV0*DTP+SMVD<*DTP2)/TP  M  85 

AMPJ=SQRT(SUMJ1**2+SUMJ2*»2)  M  86 

AMPV=SQRT{SUMV1»»2+SUMV2**2)  M  87 

AMP2  =SQRT(SUMJ12*»2+SUMJ22*»2)  M  88 

AMPV2=SQRT(SUMV12**2+SUMV22»<‘2)  M  89 

AMP=AMPJ  M  90 

PJ=ATAN2(SUMJ2,SUMJ1)*DPR  M  91 

PV=ATAN2(SUMV2,SUMV1)*DPR  M  92 

IF(PJ.LT,0.)  PJ=360.+PJ  M  93 

IF(PV.LT.0.)  PV=360,+PV  M  94 

PHASE=PJ-PV  M  95 

PJ2=ATAN2(SUMJ22,SUMJ12)*DPR  M  96 

PV2=ATAN2(SUMV22,SUMV12)*DPR  M  97 

IF(PJ2.LT.0.)  PJ2=360.+PJ2  M  50 

IF(PV2.LT.0.)  PV2=360.+PV2  M  99 

PHASE2=PJ2-PV2  M  100 

GNEG=AMP J*ABS{ COS ( PHASE/DPR ))*S/AMPV  M  101 

RFP=AMPV*»2*GNEG/2.  M  102 

DCPWR  =AMPVDC*AMPJDC*S  M  103 

EFF=RFP/DCPWR  M  104 

RETURN  M  105 

END  M  106 

BLOCK  DATA  N  1 

REAL  JPDIF,JMDIF,MUP,MUM,NOP,NOM,MUPO  ,JPO,JMAX,JDISPL,JC, lEXT, I  IN  N  2 

1T,JINT, I INTO,L,JAVG, 123, ICUR, ICURO,JMO,MJP,MJM,MJPDF,MJMDF,MJDSP,  N  3 
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2 JCSAV , JPSAV , MUMO , NP , NM , JP , JM , JMSAV , NMSAV , NPSAV , MTEM  N 

INTEGER  PRTFRQ, EXTRAP  N 

COMMON  /ABC/  MODA, MODB, MOOD, MODVP,MODVM,MODVNP,MOOVNM, MOOT, A1 ,A2,  N 
1A3,A4,A5,A6,B1,B2,B3,B4,B5,B6,D1,D2,D3,D4,D5,D6,C1,C2,C3,C4,C5,C6,  N 
2C7,C8,C9,C10,C11,C12,C13,C14.C15.C16,C17,C18,C19,C20,MUP,MUM,MUP0,  N 
3MUMO,AVP,AVM,NCP,NOM,VP1M,VM1M,CAT,CBT,PWRP,PWRM,TEMP  N 

COMMON  /CIRT/  USTAT ,U , V,DU , VO, VI , V2 ,OMEGA,OMEGA1 ,0MEGA2 ,PH 1 1 ,PH 1 2,  N 
1RD,LD,R(3),L(3) ,C( 3 ) ,CP ,T1 , I EFCY,VNEW  N 

COMMON  /DSTR/  NP(201 ) ,NM(201 ) ,DN(201 ) ,TEM(201 ) ,E{201 ) ,DE(201 ) ,VSUM  N 
1(201),VP(201),VM(201) ,JP(201) ,JM(201),JPD1F(201),JMD1F(201),JDISPL  N 
2(201 ) ,RCMBR(201 ) ,ALPHA(201 ) ,BETA(201 ) ,TPR I NT(201 )  N 

COMMON  /EFCY/  AMP, AMP2, PHASE, PJ,PV,PHASE2,PJ2,PV2, OTP, GNEG,RFP,EFF  N 
1,AMPVDC,AMPJDC,DCPWR,AMPV,AMPV2  N 

COMMON  /FUND/  M,MP1 ,N3S, PRTFRQ, NXPRNT(6) ,M0DP,M0DFCH, EXTRAP, MS, IJ,  N 
1IPOINT,FM,STEPFA,OX,DT,T,D,S,P,EL,EPSO,DIELK,DNSTY,SPCHT,THCND,DEN  N 
2SPL,DENSMI ,DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI , GAMMA, TAUSEC  N 

COMMON  /MJV/  JPO,JMO,JMAX,DJMAXM,JC,JINT,l INT,I INT0,IEXT,I23,DIDT,  N 
1DVDT,DVDT0, VC( 3 ) , VCO( 3 ) ,DVC( 3 ) ,DVCO( 3 ) , I  CUR ( 3 ) ,  I CURO( 3 ) ,MJP ,MJM,  N 
2MJPDF ,MJMDF ,MJDSP ,MTEM, VPMAX, VMMAX,CST{ 10,3 ) ,DMTEM,DD I DT ,0 1 DTO  N 

COMMON  /PLO/  TIME(401),VOLT(401),CURR(401),JAVG(401),CURR1(401),  N 
1 CURR2  ( 40 1 ) ,  CURR  3  ( 40 1 ) ,  CLTT  (401),  TAVR  (401)  .  ...N 

COMMON  /SAV/  NPSAV(201),NMSAV(201),ESAV(201),VXSAV(201),VPSAV(201)  N 

1 ,  VMSAV (201), JPSAV (201), JMSAV (201), VSAV , TSAV , JCSAV , GAMSAV , TMPSAV  N 
DATA  MODA , MODB , MODD , MOD VP , MODVM , MOD VNP , MODVNM , MOOT , A 1 , A2 , A3 , A4 , A5 ,  N 

1A6,B1,B2,B3,B4,B5,B6,01,D2,D3,D4,D5,D6,C1,C2,C3,C4,C5,C6,C7,C8,C9,  N 
2C1O,C11,C12,C13,C14,C15,C16,C17,C18,C19,C20,MUP,MUM,M'jP0,MUM0,AVP,  N 
3AVM, NOP, N0M,VP1M,VM1M, CAT, CBT,PWRP,PWRM,TEMP/8*0, 53*0.  '  N 

DATA  USTAT,U,V,DU,V0,V1,V2,0MEGA,0MEGA1,0MEGA2,PHI1,PHI2,R,L,C,CP,  N 
1T1,IEFCY/23*0.,0/  N 

DATA  AMP, AMP2, PHASE, PJ,PV,PHASE2,PJ2,PV2, DTP, GNEG,RFP,EFF,AMPVDC,  N 
1AMPJDC,DCPWR,AMPV,AMPV2/17*0./  N 

DATA  M,MP1,N3S,PRTFRQ,NXPRNT,M0DP,M0DFCH,EXTRAP,MS,!J, IPOINT,FM,ST  N 
1EPFA,DX,DT,T,D,S,P,EL,EPS0,DIELK,DNSTY,SPCHT,THCND,DENSPL,DENSMI,  N 
2DENS2,DENS2T,REC,EGAP,GAMSEC,GAMI , GAMMA, TAUSEC/ 16*0, 24*0./  N 

DATA  NP,NM,DN,TEM,E,DE,VSUM,VP,VM,JP,JM,JPDIF,JMDIF,JDISPL,RCMBR,  N 
1ALPHA,BETA,TPRINT/1918*0,/  N 

DATA  JPO , JMO , JMAX , D JMAXM , JC , J I  NT , I  I  NT , I  I NTO , I  EXT , 1 23 , D I DT , DVDT , DVD  N 
1T0,VC,VC0,DVC,DVC0, ICUR, ICURO,MJP,MJM,MJPDF,MJMDF,MJDSP, MTEM, VPMAX  N 

2, VMMAX,CST,0MTEM/70*0./  N 

DATA  TIME, VOLT, CURR, JAVG,CURR1,CURR2,CURR3/2807*0./  N 

DATA  NPSAV , NMSAV , ESAV , VXSAV , VPSAV , VMSAV , JPSAV , JMSAV , VSAV , TSAV , JCSA  N 

IV, GAMSAV, TMPSAV/813*0./  N 

END  N 
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